# Preamble ----
# File: gns_jeps_replication.R
# Description: replication for Countering Violent Etremism
# Authors: Allison Grossman, William G. Nomikos, Niloufer Siddiqui
# Last Update: 12/3/21

# **Clean workspace
rm(list=ls())

# **Working directory ----
setwd("/Users/williamnomikos/Dropbox/Work/Projects/01_Burkina Faso/JEPS/Replication")
getwd()

# **Libraries ----
#library(knitr)
#opts_chunk$set(tidy.opts=list(width.cutoff=60),tidy=TRUE)
#opts_chunk$set(echo = TRUE)
library(tidyverse)
library(magrittr)
library(dmetar) 
library(foreign)
library(kableExtra)
library(formatR)
library(cowplot)
library(TOSTER)
library(googlesheets)
library(countrycode)
library(readxl)
library(conflicted)
library(readstata13)
library(plotrix)
library(lmtest)
library(rgdal)
library(gridExtra)
library(stargazer)
library(sandwich)
library(multiwayvcov)
library(broom)
library(rstatix)
conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")

# Theme for all plots
my_theme<-theme_set(theme_bw() + theme(axis.text = element_text(size=9), plot.title = element_text(size=10, face="bold" ), axis.title=element_text(size=10, face="bold"), legend.text = element_text(size=9), legend.title = element_text(size=10, face="bold"), legend.position = "bottom"))


# ** Data Management----

  #Our main data
  data<-read_csv("gns_jeps_main.csv") #import dataset
  data_attention_fail<-data #store data frame as data_attention_fail to denote it includes those that failed attention check
  data<-subset(data,attention_fail=="No") #drop respondents who failed attention check in main data
  
  

  
  # ACLED Sample
  acled_bf<-read.csv("acled_bf.csv")
  acled_sum<-acled_bf %>% # keep only data frame of 2009-2020 of total fatalities by year
    filter(year>2009) %>% 
    group_by(year) %>% 
    summarize(total=sum(fatalities))

  # IDMC Total Displaed data
  idmc_bf<-read.csv("bf_idmc.csv")%>% 
  select(Year, Total) 

  # No displacements from year 2010-2015 so we have to add it to data as 0
  Year=c(2010:2015)
  Total=rep(0, 6)
  idmc_bf<-bind_rows(data.frame(Year,Total), idmc_bf) %>% as.data.frame()

  # Import replication sample
  data2<-read_csv("gns_jeps_sample2.csv")

#write_csv(data2,"gns_jeps_sample2.csv")

  # Create subgroups
  data_mossi<-subset(data,eth=="Mossi")
  data_nomossi<-subset(data,eth!="Mossi")
  
  data_14_plus<-subset(data,age>=14)
  
  data_muslims <- subset(data,religion=="Muslim")
  data_nonmuslims <- subset(data,religion!="Muslim")

  
  data_mossi2<-subset(data,eth=="Mossi"& religion=="Muslim")
  data_nomossi2<-subset(data,eth!="Mossi" & religion=="Muslim")
  
  # Other Samples
  
  # Afrobarometer Round 7
  r7<-read.spss("bfo_r7_data_eng.sav")
  r7<-as.data.frame(r7)
  r7_df<-r7 %>% select(Q85B, Q87A, Q87B, Q87D) # keep only these qs
  
  # Afrobarometer Round 5
  r5<-read.spss("bfo_r5_data.sav")
  r5<-as.data.frame(r5)
  r5_df<-r5 %>% select(Q78)  # Keep only question 78

  # Combine two Afrobarometer samples
  r7_df<-bind_cols(r7_df, r5_df)
  colnames(r7_df)<-c("id_strength", "nb_relig", "nb_eth", "nb_foreign", "violence_justified")
  
  # Prepare Afrobarometer Data for Comparison with our sample
  r7_df %<>% mutate(in_sample=0)
  r7_df %<>% mutate(id_strength=factor(id_strength, levels=c("I feel only (national identity)", "I feel more (national identity) than (ethnic group)", "I feel equally (national identity) and (ethnic group)", "I feel more (ethnic group) than (national identity)", "I feel only (ethnic group)"), labels=c(1:5)))
  r7_df %<>% mutate(nb_relig=factor(nb_relig, levels=c("Strongly dislike", "Somewhat dislike", "Would not care", "Somewhat like", "Strongly like"), labels=c(1:5)),
                    nb_eth=factor(nb_eth, levels=c("Strongly dislike", "Somewhat dislike", "Would not care", "Somewhat like", "Strongly like"), labels=c(1:5)), 
                    nb_foreign=factor(nb_foreign, levels=c("Strongly dislike", "Somewhat dislike", "Would not care", "Somewhat like", "Strongly like"), labels=c(1:5)))
  r7_df %<>% mutate(violence_justified=factor(violence_justified, levels=c("Agree very strongly with 1", "Agree with 1", "Agree with 2", "Agree very strongly with 2"), labels=c(1:4)))
  
  
  # Pew
  pew_df<-read.spss("pew_2015.sav")
  pew_df<-as.data.frame(pew_df)
  
  pew_df %<>% filter(COUNTRY=="Burkina Faso") # keep only BF
  pew_df %<>% select(Q23) # keep only Q23
  

  
  
# **Functions ----
#Balance Function
bal.table<-function(x,row){
  
  # Means
  mean_con <- mean(x[data$treat=="Control"], na.rm=T)
  mean_t1 <- mean(x[data$treat=="Scripture and ID"], na.rm=T)
  mean_t2 <- mean(x[data$treat=="Scripture and No ID"], na.rm=T)
  mean_t3 <- mean(x[data$treat=="Superordinate and ID"], na.rm=T)
  mean_t4 <- mean(x[data$treat=="Superordinate and No ID"], na.rm=T)
  
  # Differences
  ## Control
  diff_con_t1 <-mean_con-mean_t1
  diff_con_t2 <-mean_con-mean_t2
  diff_con_t3 <-mean_con-mean_t3
  diff_con_t4 <-mean_con-mean_t4
  ## T1
  diff_t1_t2 <-mean_t1-mean_t2
  diff_t1_t3 <-mean_t1-mean_t3
  diff_t1_t4 <-mean_t1-mean_t4
  ##T2
  diff_t2_t3 <-mean_t2-mean_t3
  diff_t2_t4 <-mean_t2-mean_t4
  ##T3
  diff_t3_t4 <-mean_t3-mean_t4
  
  
  # pvals
  ## Control
  pval_con_t1 <-as.numeric(t.test(x[data$treat=="Control"], x[data$treat=="Scripture and ID"], alternative = "two.sided", var.equal = FALSE)[3])
  pval_con_t2 <-as.numeric(t.test(x[data$treat=="Control"], x[data$treat=="Scripture and No ID"], alternative = "two.sided", var.equal = FALSE)[3])
  pval_con_t3 <-as.numeric(t.test(x[data$treat=="Control"], x[data$treat=="Superordinate and ID"], alternative = "two.sided", var.equal = FALSE)[3])
  pval_con_t4 <-as.numeric(t.test(x[data$treat=="Control"], x[data$treat=="Superordinate and No ID"], alternative = "two.sided", var.equal = FALSE)[3])
  ## T1
  pval_t1_t2 <-as.numeric(t.test(x[data$treat=="Scripture and ID"], x[data$treat=="Scripture and No ID"], alternative = "two.sided", var.equal = FALSE)[3])
  pval_t1_t3 <-as.numeric(t.test(x[data$treat=="Scripture and ID"], x[data$treat=="Superordinate and ID"], alternative = "two.sided", var.equal = FALSE)[3])
  pval_t1_t4 <-as.numeric(t.test(x[data$treat=="Scripture and ID"], x[data$treat=="Superordinate and No ID"], alternative = "two.sided", var.equal = FALSE)[3])
  ##T2
  pval_t2_t3 <-as.numeric(t.test(x[data$treat=="Scripture and No ID"], x[data$treat=="Superordinate and ID"], alternative = "two.sided", var.equal = FALSE)[3])
  pval_t2_t4 <-as.numeric(t.test(x[data$treat=="Scripture and No ID"], x[data$treat=="Superordinate and No ID"], alternative = "two.sided", var.equal = FALSE)[3])
  ##T3
  pval_t3_t4 <-as.numeric(t.test(x[data$treat=="Superordinate and ID"], x[data$treat=="Superordinate and No ID"], alternative = "two.sided", var.equal = FALSE)[3])
  
  
  t <- round(cbind(
    mean_con,mean_t1,mean_t2,mean_t3,mean_t4,
    pval_con_t1,pval_con_t2,pval_con_t3,pval_con_t4,
    pval_t1_t2,pval_t1_t3,pval_t1_t4,
    pval_t2_t3,pval_t2_t4,
    pval_t3_t4
  ),digits=2)
  tab_balance[row,]<-t
  tab_balance<<-tab_balance
}
# Balance Function for conparing different samples
bal.table_sample<-function(x,row){
  
  # Means
  mean_sample1 <- mean(x[summary_data_all$sample=="Main"], na.rm=T)
  mean_sample2 <- mean(x[summary_data_all$sample=="Replication"], na.rm=T)
  
  
  # Differences
  ## Control
  diff <-mean_sample1-mean_sample2
  
  # pvals
  ## Control
  pval_con_t1 <-as.numeric(t.test(x[summary_data_all$sample=="Main"], x[summary_data_all$sample=="Replication"], alternative = "two.sided", var.equal = FALSE)[3])
  
  
  t <- round(cbind(
    mean_sample1, mean_sample2, diff, pval_con_t1),digits=2)
  tab_balance[row,]<-t
  tab_balance<<-tab_balance
}


# Results in Text ----
# These are the results we discuss in the body of the manuscript
  # Treatment effect on tolerance of members of other ethnic groups (full sample)
  summary(lm(ethnic_binary ~ treat_binary, data=data)) # regression results
  confint(lm(ethnic_binary ~ treat_binary, data=data))[2,] # 95% CI only
  
  # Treatment effect on trust of members of other ethnic groups (full sample)
  summary(lm(trust_eth_binary ~ treat_binary, data=data))  # regression results
  confint(lm(trust_eth_binary ~ treat_binary, data=data))[2,]  # 95% CI only

  # Treatment effect on tolerance of members of other ethnic groups (Mossi sample)
  summary(lm(ethnic_binary ~ treat_binary, data=data_mossi)) # regression results
  confint(lm(ethnic_binary ~ treat_binary, data=data_mossi))[2,] # 95% CI only

  # Treatment effect on trust of members of other ethnic groups (Mossi sample)
  summary(lm(trust_eth_binary ~ treat_binary, data=data_mossi))
  confint(lm(trust_eth_binary ~ treat_binary, data=data_mossi))[2,]
  
  # Treatment effect on trust of members of other religious groups (Mossi sample)
  summary(lm(trust_rel_binary ~ treat_binary, data=data_mossi))
  confint(lm(trust_rel_binary ~ treat_binary, data=data_mossi))[2,]
  
  # Treatment effect on trust of members of other religious groups (non-Mossi sample)
  summary(lm(trust_rel_binary ~ treat_binary, data=data_nomossi))
  confint(lm(trust_rel_binary ~ treat_binary, data=data_nomossi))[2,]

  # Treatment effect on willingness to send message to government (Mossi sample)
  summary(lm(msg_gov ~ treat_binary, data=data_mossi))
  confint(lm(msg_gov ~ treat_binary, data=data_mossi))[2,]


  # Treatment effect on willingness to send message to government (non-Mossi sample)
  summary(lm(msg_gov ~ treat_binary, data=data_nomossi))
  confint(lm(msg_gov ~ treat_binary, data=data_nomossi))[2,]

# Figures ----

  # ** Figure 1 ----
    # First panel: Fatalities
    acled_plot<-acled_sum %>% ggplot(aes(x=year, y=total)) + geom_line(color="red") + geom_vline(xintercept=2016, linetype=2) +
    theme_bw() + theme(axis.text = element_text(size=9), plot.title = element_text(size=12, face="bold"), axis.title=element_text(size=10, face="bold"), legend.text = element_text(size=9), legend.title = element_text(size=10, face="bold"), legend.position = "bottom") + scale_x_continuous(breaks=c(2010, 2012, 2014, 2016, 2018, 2020)) + scale_y_continuous(limits=c(0, 3000), breaks=c(seq(0, 3000, by=500))) + labs(y="Fatalities from violent events (ACLED)", x="Year")
  
    # Second panel: Displacement
    idmc_plot<-idmc_bf %>% ggplot(aes(x=Year, y=Total)) + geom_line(color="orange") + geom_vline(xintercept=2016, linetype=2) +
    theme_bw() + theme(axis.text = element_text(size=9), plot.title = element_text(size=12, face="bold"), axis.title=element_text(size=10, face="bold"), legend.text = element_text(size=9), legend.title = element_text(size=10, face="bold"), legend.position = "bottom") + scale_x_continuous(breaks=c(2010, 2012, 2014, 2016, 2018, 2020)) + scale_y_continuous(limits=c(0, 600000), breaks=c(seq(0, 600000, by=100000))) + labs(y="Internally displaced persons (IDMC)", x="Year")
  
    #Combine two panels
    fig1<-grid.arrange(acled_plot, idmc_plot, ncol=2)
    
    #Save
    save_plot("fig1.png", fig1, base_width = 6, base_height = 4)
  
  # ** Figure 2 ----
    # (a) Attitudes toward outgroups
    # Regress each outcome on the binary treatment indicator and store the results
    fit1<-lm(religion_binary ~ treat_binary, data=data)
    fit2<-lm(ethnic_binary ~ treat_binary, data=data)
    fit3<-lm(trust_rel_binary ~ treat_binary, data=data)
    fit4<-lm(trust_eth_binary ~ treat_binary, data=data)
    reg_list<-list(fit1, fit2, fit3, fit4)
        
    # Grab the coefficients, 95%, and 90% CI of each and store in a dataframe
    fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
    fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
    fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3), confint(fit3, level=.90))
    fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4), confint(fit4, level=.90))
    
    # Create columns that label treatment and outcomes 
    fit1.df$names<-rownames(fit1.df)
    fit1.df$outcome<-"Tolerance of neighbor\nfrom  different religion"
    fit2.df$names<-rownames(fit2.df)
    fit2.df$outcome<-"Tolerance of neighbor\nfrom  different ethnicity"
    fit3.df$names<-rownames(fit3.df)
    fit3.df$outcome<-"Trusts members of \nother religious groups"
    fit4.df$names<-rownames(fit4.df)
    fit4.df$outcome<-"Trusts members of \nother ethnic groups"
    
    #Create labels for the remaining columns of dataframe
    colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit3.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit4.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    
    #Bring together the four regressions
    all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df)
    all_reg %<>% filter(!names=="(Intercept)") #keep only treatment effect
    all_reg$order <- seq(from = nrow(all_reg), to = 1, by = -1) # create a counter to reverse order of presentation
    
    #Create figure and store it
    fig2a<-ggplot(all_reg,aes(x=reorder(outcome,order),y=coef))+
      coord_flip()+
      geom_hline(yintercept = 0, color = "gray80", lty = "dashed")+
      geom_point(position=position_dodge(width = .5), size = 2.5) +
      theme(panel.background = element_blank(),
            legend.title = element_blank(), 
            plot.title = element_text(size = 12),
            panel.border = element_rect(colour = "gray", fill=NA, size=.8),
            legend.position = "bottom",
            axis.title.x  = element_text(size=10),
            axis.text.x  = element_text(angle=0, vjust=.5, hjust = .5, size=8))+
      geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
                    position=position_dodge(width = 0.5),
                    width=0, lwd = 0.4)+
      geom_errorbar(aes(ymin=lower_90, ymax=upper_90), 
                    position=position_dodge(width = 0.5),
                    width=0, lwd = 1) +
      scale_y_continuous(limits=c(-0.15,0.15)) + 
      labs(x="", y="Average Treatment Effect",title="(a) Attitudes toward Outgroups")

    # (b) Attitudes about extremism
    # Repeat same steps as above for attitudes about extremism
    fit1<-lm(concern_binary ~ treat_binary, data=data)
    fit2<-lm(justified_binary ~ treat_binary, data=data)
    
    reg_list<-list(fit1, fit2)
    
    fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
    fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
    
    
    fit1.df$names<-rownames(fit1.df)
    fit1.df$outcome<-"Concerned about\nviolent extremism"
    fit2.df$names<-rownames(fit2.df)
    fit2.df$outcome<-"Believes violence\n can be justified"
    
    colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    
    all_reg<-bind_rows(fit1.df, fit2.df)
    all_reg %<>% filter(!names=="(Intercept)") 
    all_reg$order <- seq(from = nrow(all_reg), to = 1, by = -1)
    
    fig2b<-ggplot(all_reg,aes(x=reorder(outcome,order),y=coef))+
      coord_flip()+
      geom_hline(yintercept = 0, color = "gray80", lty = "dashed")+
      geom_point(position=position_dodge(width = .5), size = 2.5) +
      theme(panel.background = element_blank(),
            legend.title = element_blank(), 
            plot.title = element_text(size = 10),
            panel.border = element_rect(colour = "gray", fill=NA, size=.8),
            legend.position = "bottom",
            axis.title.x  = element_text(size=8),
            axis.text.x  = element_text(angle=0, vjust=.5, hjust = .5, size=6))+
      geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
                    position=position_dodge(width = 0.5),
                    width=0, lwd = 0.4)+
      geom_errorbar(aes(ymin=lower_90, ymax=upper_90), 
                    position=position_dodge(width = 0.5),
                    width=0, lwd = 1) +
      scale_y_continuous(limits=c(-0.15,0.15)) + 
      labs(x="", y="Average Treatment Effect",title="(b) Attitudes about Extremism")

    


    # (c) Civic Participation
    # Repeat same steps as above for civic participation
        
    fit1<-lm(message_binary ~ treat_binary, data=data) #any message
    fit2<-lm(msg_gov~ treat_binary, data=data) # Message to government
    fit3<-lm(msg_islamic ~ treat_binary, data=data) # Message to Islamic council
    
    
    reg_list<-list(fit1, fit2, fit3)
    
    fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
    fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
    fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3), confint(fit3, level=.90))
    
    fit1.df$names<-rownames(fit1.df)
    fit1.df$outcome<-"Sends any message"
    
    fit2.df$names<-rownames(fit2.df)
    fit2.df$outcome<-"Sends message to\ngovernment"
    
    fit3.df$names<-rownames(fit3.df)
    fit3.df$outcome<-"Sends message to\nIslamic council"
    
    
    colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit3.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    
    
    all_reg<-bind_rows(fit1.df, fit2.df, fit3.df)
    all_reg %<>% filter(!names=="(Intercept)") 
    all_reg$order <- seq(from = nrow(all_reg), to = 1, by = -1)
    
    
    fig2c<-ggplot(all_reg,aes(x=reorder(outcome,order),y=coef))+
      coord_flip()+
      geom_hline(yintercept = 0, color = "gray80", lty = "dashed")+
      geom_point(position=position_dodge(width = .5), size = 2.5) +
      theme(panel.background = element_blank(),
            legend.title = element_blank(), 
            plot.title = element_text(size = 10),
            panel.border = element_rect(colour = "gray", fill=NA, size=.8),
            legend.position = "bottom",
            axis.title.x  = element_text(size=8),
            axis.text.x  = element_text(angle=0, vjust=.5, hjust = .5, size=6))+
      geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
                    position=position_dodge(width = 0.5),
                    width=0, lwd = 0.4)+
      geom_errorbar(aes(ymin=lower_90, ymax=upper_90), 
                    position=position_dodge(width = 0.5),
                    width=0, lwd = 1) +
      scale_y_continuous(limits=c(-0.15,0.15)) + 
      labs(x="", y="Average Treatment Effect",title="(c) Civic Participation")



    # Combine figures and save
    fig2<-grid.arrange(fig2a, arrangeGrob(fig2b, fig2c), ncol = 2,widths=c(40,35))
    ggsave("fig2.png", fig2, width = 8, height = 4)



  # ** Figure 3 ----
  # Figure 3 is created using the same procedured as Figure except using subgroups of Mossi/non-Mossi
    # (a) Attitudes toward outgroups
    fit1<-lm(religion_binary ~ treat_binary, data=data_mossi)
    fit2<-lm(religion_binary ~ treat_binary, data=data_nomossi)
    fit3<-lm(ethnic_binary ~ treat_binary, data=data_mossi)
    fit4<-lm(ethnic_binary ~ treat_binary, data=data_nomossi)
    fit5<-lm(trust_rel_binary ~ treat_binary, data=data_mossi)
    fit6<-lm(trust_rel_binary ~ treat_binary, data=data_nomossi)
    fit7<-lm(trust_eth_binary ~ treat_binary, data=data_mossi)
    fit8<-lm(trust_eth_binary ~ treat_binary, data=data_nomossi)
    
    reg_list<-list(fit1, fit2, fit3, fit4,fit5,fit6,fit7,fit8)
    
    fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
    fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
    fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3), confint(fit3, level=.90))
    fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4), confint(fit4, level=.90))
    fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5), confint(fit5, level=.90))
    fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6), confint(fit6, level=.90))
    fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7), confint(fit7, level=.90))
    fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8), confint(fit8, level=.90))
    
    fit1.df$names<-rownames(fit1.df)
    fit1.df$outcome<-"Tolerance of neighbor\nfrom  different religion"
    fit2.df$names<-rownames(fit2.df)
    fit2.df$outcome<-"Tolerance of neighbor\nfrom  different religion"
    fit3.df$names<-rownames(fit3.df)
    fit3.df$outcome<-"Tolerance of neighbor\nfrom  different ethnicity"
    fit4.df$names<-rownames(fit4.df)
    fit4.df$outcome<-"Tolerance of neighbor\nfrom  different ethnicity"
    fit5.df$names<-rownames(fit5.df)
    fit5.df$outcome<-"Trusts members of \nother religious groups"
    fit6.df$names<-rownames(fit6.df)
    fit6.df$outcome<-"Trusts members of \nother religious groups"
    fit7.df$names<-rownames(fit7.df)
    fit7.df$outcome<-"Trusts members of \nother ethnic groups"
    fit8.df$names<-rownames(fit8.df)
    fit8.df$outcome<-"Trusts members of \nother ethnic groups"
    
    
    
    colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit3.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit4.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    
    colnames(fit5.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit6.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit7.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit8.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    
    all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df,fit5.df,fit6.df,fit7.df,fit8.df)
    all_reg %<>% filter(!names=="(Intercept)") 
    all_reg$names<-c("Mossi","Non-Mossi","Mossi","Non-Mossi","Mossi","Non-Mossi","Mossi","Non-Mossi")
    all_reg$order <- seq(from = nrow(all_reg), to = 1, by = -1)
    
    
    fig3a<-ggplot(all_reg,aes(x=reorder(outcome,order),y=coef,group = names,color = names,shape=names))+
      coord_flip()+
      geom_hline(yintercept = 0, color = "gray80", lty = "dashed")+
      geom_point(position=position_dodge(width = .5), size = 2.5) +
      theme(panel.background = element_blank(),
            legend.title = element_blank(), 
            plot.title = element_text(size = 10),
            panel.border = element_rect(colour = "gray", fill=NA, size=.8),
            legend.position = "bottom",
            axis.title.x  = element_text(size=9),
            axis.text.x  = element_text(angle=0, vjust=.5, hjust = .5, size=8))+
      scale_colour_manual(values=c("#ba0000", "#0072B2"))+
      geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
                    position=position_dodge(width = 0.5),
                    width=0, lwd = 0.4)+
      geom_errorbar(aes(ymin=lower_90, ymax=upper_90), 
                    position=position_dodge(width = 0.5),
                    width=0, lwd = 1) +
      scale_y_continuous(limits=c(-0.4,0.4)) + 
      labs(x="", y="Average Treatment Effect",title="(a) Attitudes toward Outgroups")

    
    
    #(b) Attitudes about extremism
    fit1<-lm(concern_binary ~ treat_binary, data=data_mossi)
    fit2<-lm(concern_binary ~ treat_binary, data=data_nomossi)
    fit3<-lm(justified_binary ~ treat_binary, data=data_mossi)
    fit4<-lm(justified_binary ~ treat_binary, data=data_nomossi)
    
    
    reg_list<-list(fit1, fit2, fit3, fit4)
    
    fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
    fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
    fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3), confint(fit3, level=.90))
    fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4), confint(fit4, level=.90))
    
    
    fit1.df$names<-rownames(fit1.df)
    fit1.df$outcome<-"Concerned about \nviolent extremism"
    fit2.df$names<-rownames(fit2.df)
    fit2.df$outcome<-"Concerned about \nviolent extremism"
    fit3.df$names<-rownames(fit3.df)
    fit3.df$outcome<-"Believes violence\ncan be justified"
    fit4.df$names<-rownames(fit4.df)
    fit4.df$outcome<-"Believes violence\ncan be justified"
    
    colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit3.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit4.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    
    all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df)
    all_reg %<>% filter(!names=="(Intercept)") 
    all_reg$names<-c("Mossi","Non-Mossi","Mossi","Non-Mossi")
    all_reg$order <- seq(from = nrow(all_reg), to = 1, by = -1)
    
    fig3b<-ggplot(all_reg,aes(x=reorder(outcome,order),y=coef,group = names,color = names,shape=names))+
      coord_flip()+
      geom_hline(yintercept = 0, color = "gray80", lty = "dashed")+
      geom_point(position=position_dodge(width = .5), size = 2.5) +
      theme(panel.background = element_blank(),
            legend.title = element_blank(), 
            plot.title = element_text(size = 10),
            panel.border = element_rect(colour = "gray", fill=NA, size=.8),
            legend.position = "bottom",
            axis.title.x  = element_text(size=8),
            axis.text.y  = element_text(size=7),
            axis.text.x  = element_text(angle=0, vjust=.5, hjust = .5, size=6))+
      scale_colour_manual(values=c("#ba0000", "#0072B2"))+
      geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
                    position=position_dodge(width = 0.5),
                    width=0, lwd = 0.4)+
      geom_errorbar(aes(ymin=lower_90, ymax=upper_90), 
                    position=position_dodge(width = 0.5),
                    width=0, lwd = 1) +
      scale_y_continuous(limits=c(-0.4,0.4)) + 
      labs(x="", y="Average Treatment Effect",title="(b) Attitudes about Extremism")


    # (c) Civic Participation
    fit1<-lm(message_binary ~ treat_binary, data=data_mossi)
    fit2<-lm(message_binary ~ treat_binary, data=data_nomossi)
    fit3<-lm(msg_gov ~ treat_binary, data=data_mossi)
    fit4<-lm(msg_gov~ treat_binary, data=data_nomossi)
    fit5<-lm(msg_islamic ~ treat_binary, data=data_mossi)
    fit6<-lm(msg_islamic~ treat_binary, data=data_nomossi)
    
    
    reg_list<-list(fit1, fit2, fit3, fit4,fit5,fit6)
    
    fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
    fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
    fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3), confint(fit3, level=.90))
    fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4), confint(fit4, level=.90)) 
    fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5), confint(fit5, level=.90))
    fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6), confint(fit6, level=.90))
    
    
    fit1.df$names<-rownames(fit1.df)
    fit1.df$outcome<-"Sends any message"
    fit2.df$names<-rownames(fit2.df)
    fit2.df$outcome<-"Sends any message"
    fit3.df$names<-rownames(fit3.df)
    fit3.df$outcome<-"Sends message to\ngovernment"
    fit4.df$names<-rownames(fit4.df)
    fit4.df$outcome<-"Sends message to\ngovernment"
    fit5.df$names<-rownames(fit5.df)
    fit5.df$outcome<-"Sends message to\nIslamic council"
    fit6.df$names<-rownames(fit6.df)
    fit6.df$outcome<-"Sends message to\nIslamic council"
    
    colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit3.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit4.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    
    colnames(fit5.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    colnames(fit6.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
    
    
    all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df,fit5.df,fit6.df)
    all_reg %<>% filter(!names=="(Intercept)") 
    all_reg$names<-c("Mossi","Non-Mossi","Mossi","Non-Mossi","Mossi","Non-Mossi")
    all_reg$order <- seq(from = nrow(all_reg), to = 1, by = -1)
    
    
    fig3c<-ggplot(all_reg,aes(x=reorder(outcome,order),y=coef,group = names,color = names,shape=names))+
      coord_flip()+
      geom_hline(yintercept = 0, color = "gray80", lty = "dashed")+
      geom_point(position=position_dodge(width = .5), size = 2.5) +
      theme(panel.background = element_blank(),
            legend.title = element_blank(), 
            plot.title = element_text(size = 10),
            panel.border = element_rect(colour = "gray", fill=NA, size=.8),
            legend.position = "bottom",
            axis.title.x  = element_text(size=8),
            axis.text.x  = element_text(angle=0, vjust=.5, hjust = .5, size=6))+
      scale_colour_manual(values=c("#ba0000", "#0072B2"))+
      geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
                    position=position_dodge(width = 0.5),
                    width=0, lwd = 0.4)+
      geom_errorbar(aes(ymin=lower_90, ymax=upper_90), 
                    position=position_dodge(width = 0.5),
                    width=0, lwd = 1) +
      scale_y_continuous(limits=c(-0.4,0.4)) + 
      labs(x="", y="Average Treatment Effect",title="(c) Civic Participation")
    
  # Combine figures, store, and save
  fig3<-grid.arrange(fig3a, arrangeGrob(fig3b, fig3c), ncol = 2,widths=c(40,35))
  ggsave("fig3.png", fig3, width = 8, height = 4)
  
  





  
  
  
  
  
# Appendix Tables ----
  

  names(data)
  # ** Table C1 ----
  # This is a summary statistics table for the main study
  # Select only the variables to be summarized from the data
  summary_data<-data %>% 
    select(age, Muslim, Ouaga, Girl, econ_hh ,laws_qran, pray_often, 
           relig_important, violence_exp, violence_heard,join_econ, diff_eth, diff_relig) %>% 
    as.data.frame 
  # We use the stargazer package to produce the Table directly in LaTeX below:
  summary_data %>% 
    stargazer(no.space = T,
              summary.stat = c("n", "mean", "sd", "median", "min", "max"), 
              covariate.labels = c("Age", "Muslim", "Live in Ouagadougou", "Girl", "Household economic status", 
                                   "Laws should strictly follow Qu'ran", "Pray more than once a week", "Religion very important",
                                   "Experienced violence in past year", "Heard of violence in past year", 
                                   "People join extremists because lack of job/money", "Relations between different ethnicities", 
                                   "Relations between different religions") )

  # ** Table C2 ----
  # This is a summary statistics table for the replication study
  # Select only the variables to be summarized from the data
  summary_data<-data2 %>% 
    select(age, Muslim, Girl, econ_hh ,laws_qran, pray_often, 
           relig_important, violence_exp, violence_heard,join_econ, diff_eth, diff_relig) %>% 
    as.data.frame 
  # We use the stargazer package to produce the Table directly in LaTeX below:
  summary_data %>% 
    stargazer(no.space = T,
              summary.stat = c("n", "mean", "sd", "median", "min", "max"), 
              covariate.labels = c("Age", "Muslim", "Girl", "Household economic status", 
                                   "Laws should strictly follow Qu'ran", "Pray more than once a week", "Religion very important",
                                   "Experienced violence in past year", "Heard of violence in past year", 
                                   "People join extremists because lack of job/money", "Relations between different ethnicities", 
                                   "Relations between different religions") )
  
  
  # ** Table C3 ----
  # Balance on Demographic Covariates between Treatments (means)
  
  # Make sure that you've run bal.table, a function we defined
  # First create an empty table, which stores all of the means and pvals the function produces
  # We have 13 rows for 13 variables of interest and 15 columns for means and pvals
  tab_balance<-matrix(NA,13,15)
  
  # Next, run the function for each variables of interest, indicating which row of the table to store it
  bal.table(data$age,1)
  bal.table(data$Muslim,2)
  bal.table(data$Ouaga,3)
  bal.table(data$Girl,4)
  bal.table(data$econ_hh,5)
  bal.table(data$laws_qran,6)
  bal.table(data$pray_often,7)
  bal.table(data$relig_important,8)
  bal.table(data$violence_exp,9)
  bal.table(data$violence_heard,10)
  bal.table(data$join_econ,11)
  bal.table(data$diff_eth,12)
  bal.table(data$diff_relig,13)
  
  #Name rows after each variable
  rownames(tab_balance)<-c("Age", "Muslim", "Ouagadougou", "Girl", "Household econ.", "Laws should follow Qu'ran", "Pray frequent", "Religion  important","Experienced violence ", "Heard of violence ", "Join extremists bc money", "Relations  diff eth", "Relations  diff rel")
  
  #Label columns
  colnames(tab_balance) <- c(
    "mean_con","mean_t1","mean_t2","mean_t3","mean_t4",
    "pval_con_t1","pval_con_t2","pval_con_t3","pval_con_t4",
    "pval_t1_t2","pval_t1_t3","pval_t1_t4",
    "pval_t2_t3","pval_t2_t4",
    "pval_t3_t4"
  )
  
  # Create a table that stores the first five columns (the means)
  tab_balance_means<-tab_balance[,1:5]
  # Print the table in R
  tab_balance_means
  # Create the LaTeX code
  kable(tab_balance_means,"latex",booktabs=T,linesep = "",caption = "Balance on Demographic Covariates between Treatments (means)") %>%
    kable_styling(latex_options = "striped", font_size = 10)
  
  
  # ** Table C4 ----
  # Balance on Demographic Covariates between Treatments (pvals)
  # This assumes you already ran all the code for Table C3
  # Create a table that stores the last nine columns of the table(the pvalues)
  tab_balance_pvals<-tab_balance[,6:15]
  # Print in R
  tab_balance_pvals
  # Create the LaTeX code
  kable(tab_balance_2,"latex",booktabs=T,linesep = "",caption = "Balance on Demographic Covariates between Treatments (pvals)") %>%
    kable_styling(latex_options = "striped", font_size = 8)
  
  # ** Table C5 ----
  # Balance on demographic covariates between main and replication sample students in Kenedougou
  # Combine main and replication data
  data_all<-bind_rows(data, data2)

  #subset to only Kenedougou respondents 
  data_all %<>% filter(location=="KENEDOUGOU") 

  
  table(summary_data_all$sample)
  # Put all covariates of interest in one data frame
  summary_data_all<-data_all %>% 
    select(age, Muslim, Girl, econ_hh ,laws_qran, pray_often, relig_important, 
           violence_exp, violence_heard,join_econ, diff_eth, diff_relig, sample) %>% 
    as.data.frame 

  # Make sure that you've run bal.table_sample, a function we defined
  # First create an empty table, which stores all of the means and pvals the function produces
  # We have 13 rows for 12 variables of interest and 4 columns for means, diff, and pvals
  tab_balance<-matrix(NA,12,4)
  # Next, run the function for each variables of interest, indicating which row of the table to store it
  bal.table_sample(summary_data_all$age,1)
  bal.table_sample(summary_data_all$age,1)
  bal.table_sample(summary_data_all$Muslim,2)
  bal.table_sample(summary_data_all$Girl,3)
  bal.table_sample(summary_data_all$econ_hh,4)
  bal.table_sample(summary_data_all$laws_qran,5)
  bal.table_sample(summary_data_all$pray_often,6)
  bal.table_sample(summary_data_all$relig_important,7)
  bal.table_sample(summary_data_all$violence_exp,8)
  bal.table_sample(summary_data_all$violence_heard,9)
  bal.table_sample(summary_data_all$join_econ,10)
  bal.table_sample(summary_data_all$diff_eth,11)
  bal.table_sample(summary_data_all$diff_relig,12)
  
  # Row Labels
  rownames(tab_balance)<-c("Age", "Muslim", "Girl", "Household econ.", "Laws should follow Qu'ran", "Pray frequent", "Religion  important","Experienced violence ", "Heard of violence ", "Join extremists bc money", "Relations diff eth", "Relations diff rel")
  # Column labels
  colnames(tab_balance)<-c("Experimental sample", "Replication sample", "Difference", "p-value")
  # Print in R
  tab_balance
  
  # Print in LaTeX
  kable(tab_balance,"latex",booktabs=T,linesep = "", # latex settings
        caption = "Balance on Demographic Covariates between Kenedougou samples" # caption
  ) %>%
    kable_styling(latex_options = "striped", font_size = 10)
  
# ** Table C6 ----
  
  #Pull questions that are in afrobarometer
  sample_df<-data %>% 
    select("id_strength", "nb_relig", "nb_eth", "nb_foreign", "violence_justified")
  # Create indicator that it's in or sample 
  sample_df %<>% mutate(in_sample=1)
  
  # Make variables from our questionnaire into factor variables
  sample_df %<>% mutate(id_strength=factor(id_strength, levels=c("I feel only Burkinabe", "I feel more Burkinabe than [respondent’s ethnic group]", "I feel equally Burkinabe and [respondent’s ethnic group]","I feel more [respondent’s ethnic group than Burkinabe", "I feel only [respondent’s ethnic group]"), labels=c(1:5)))
  sample_df %<>% mutate(nb_relig=factor(nb_relig, levels=c("Strongly dislike", "Somewhat dislike", "Neither like nor dislike", "Somewhat like", "Strongly like"), labels=c(1:5)),
                        nb_eth=factor(nb_eth, levels=c("Strongly dislike", "Somewhat dislike", "Neither like nor dislike", "Somewhat like", "Strongly like"), labels=c(1:5)), 
                        nb_foreign=factor(nb_foreign, levels=c("Strongly dislike", "Somewhat dislike", "Neither like nor dislike", "Somewhat like", "Strongly like"), labels=c(1:5)))
  sample_df %<>% mutate(violence_justified=factor(violence_justified, levels=c(1:4), labels=c(1:4)))
  
  #combine data sets
  temp_df<-bind_rows(r7_df, sample_df)
  temp_df %<>% mutate(id_strength=as.numeric(as.character(id_strength)), 
                      nb_relig=as.numeric(as.character(nb_relig)), 
                      nb_eth=as.numeric(as.character(nb_eth)),
                      nb_foreign=as.numeric(as.character(nb_foreign)), 
                      violence_justified=as.numeric(as.character(violence_justified)))
  
  # Create data frame that includes results of ttests
  t_df<-temp_df %>% 
    pivot_longer(-in_sample, names_to = "variables", values_to = "value") %>% 
    group_by(variables) %>%
    t_test(value ~ in_sample, detailed=T) %>% 
    as.data.frame()
  
  # Select list of variables in the comparison and print table in Latex
  t_df %<>% select(variables, estimate, estimate1, estimate2, p) %>% mutate(across(estimate:p, ~ round(., digits=3)))
  rm(temp_df)
  colnames(t_df)<-c("variable", "difference", "afrobarometer", "sample", "p_value")
  t_df$variable<-c("Strength of  ethnic  identity", "Trust neighbor of different ethnicity", "Trust foreign neighbor", "Trust neighbor of different religion", "Use of violence is justified")
  t_df %>% kable("latex",booktabs=T,linesep = "",caption = "T-tests comparing sample means for Afrobarometer and experimetnal samples")

  
  # ** Table C7----
  # Pull question that is in our questionnaire and Pew
  sample_df2<-data %>% select(concern)
  # Binary indicator for sample
  sample_df2 %<>% mutate(in_sample=1)
  # Reshape pew data frame
  pew_df %<>% mutate(concern=factor(Q23, levels=c("Not at all concerned",  "Not too concerned",  "Somewhat concerned", "Very concerned"), labels=c(1:4))) %>% select(concern)  %>% mutate(concern=as.numeric(as.character(concern))) %>% mutate(in_sample=0)
  # Combine samples
  temp_df2<-bind_rows(sample_df2, pew_df)
  
  # Create data frame that includes results of ttests
  t_df2<-temp_df2 %>% pivot_longer(-in_sample, names_to = "variables", values_to = "value") %>% group_by(variables) %>%
    t_test(value ~ in_sample, detailed=T) %>% as.data.frame() %>% select(variables, estimate, estimate1, estimate2, p) %>% mutate(across(estimate:p, ~ round(., digits=3)))
  colnames(t_df2)<-c("variable", "difference", "pew", "sample", "p_value")
  
  # Label variable
  t_df2$variable<-c("Concern about spread of extremism")
  # Create label
  t_df2 %>% kable("latex",booktabs=T,linesep = "",caption = "T-tests comparing sample means for Pew (2015) and experimental samples")
  
  
  
# Appendix Robustness Checks ----
  
  # ** Permutation Tests ----
  
  # Set 2x2 panel
  par(mfrow=c(2,2))
  
  set.seed(06511) #for reproducibility
  nsim <- 1000
  res <- numeric(nsim) # set aside space for results
  
  # Religion Binary
  for (i in 1:nsim) {
    # standard approach: scramble response value
    perm <- sample(nrow(data))
    bdat <- transform(data,treat_binary=treat_binary[perm])
    # compute & store difference in means; store the value
    res[i] <- mean(bdat[bdat$treat_binary==1,"religion_binary"])-
      mean(bdat[bdat$treat_binary==0,"religion_binary"])
  }
  obs <- mean(data$religion_binary[data$treat_binary==1])-
    mean(data$religion_binary[data$treat_binary==0])
  res <- c(res,obs)
  ## Plot
  hist(res,col="white",main="Neighbor of Different Religion",cex.main=1,cex.lab=1,freq=T,breaks=30,xlim=c(-0.10,0.10),ylim=c(0,200),xlab="Difference of Means")
  abline(v=obs,col="red")
  pval<-round(mean(abs(res)>=abs(obs)),3)
  text(x=0.05,y=150,label=paste("D=",round(obs,3),"\n p=",pval,sep = ""),cex=0.8)
  
  # Ethnic Binary
  for (i in 1:nsim) {
    ## standard approach: scramble response value
    perm <- sample(nrow(data))
    bdat <- transform(data,treat_binary=treat_binary[perm])
    ## compute & store difference in means; store the value
    res[i] <- mean(bdat[bdat$treat_binary==1,"ethnic_binary"])-
      mean(bdat[bdat$treat_binary==0,"ethnic_binary"])
  }
  obs <- mean(data$ethnic_binary[data$treat_binary==1])-
    mean(data$ethnic_binary[data$treat_binary==0])
  res <- c(res,obs)
  ## Plot
  hist(res,col="white",main="Neighbor of Different Ethnicity",cex.main=1,cex.lab=1, freq=T,breaks=30,xlim=c(-0.10,0.10),ylim=c(0,200),xlab="Difference of Means")
  abline(v=obs,col="red")
  pval<-round(mean(abs(res)>=abs(obs)),3)
  text(x=0.05,y=150,label=paste("D=",round(obs,3),"\n p=",pval,sep = ""),cex=0.8)
  
  # Trust Religio Other
  for (i in 1:nsim) {
    ## standard approach: scramble response value
    perm <- sample(nrow(data))
    bdat <- transform(data,treat_binary=treat_binary[perm])
    ## compute & store difference in means; store the value
    res[i] <- mean(bdat[bdat$treat_binary==1,"trust_rel_binary"])-
      mean(bdat[bdat$treat_binary==0,"trust_rel_binary"])
  }
  obs <- mean(data$trust_rel_binary[data$treat_binary==1])-
    mean(data$trust_rel_binary[data$treat_binary==0])
  res <- c(res,obs)
  
  # Plot
  hist(res,col="white",main="Trust Other Religious Group",cex.main=1,cex.lab=1, freq=T,breaks=30,xlim=c(-0.10,0.10),ylim=c(0,200),xlab="Difference of Means")
  abline(v=obs,col="red")
  pval<-round(mean(abs(res)>=abs(obs)),3)
  text(x=0.05,y=150,label=paste("D=",round(obs,3),"\n p=",pval,sep = ""),cex=0.8)
  
  
  # Ethnic Binary
  for (i in 1:nsim) {
    ## standard approach: scramble response value
    perm <- sample(nrow(data))
    bdat <- transform(data,treat_binary=treat_binary[perm])
    ## compute & store difference in means; store the value
    res[i] <- mean(bdat[bdat$treat_binary==1,"trust_eth_binary"])-
      mean(bdat[bdat$treat_binary==0,"trust_eth_binary"])
  }
  obs <- mean(data$trust_eth_binary[data$treat_binary==1])-
    mean(data$trust_eth_binary[data$treat_binary==0])
  res <- c(res,obs)
  
  # Plot
  hist(res,col="white",main="Trust Other Ethnic Group",cex.main=1,
       cex.lab=1, freq=T,breaks=30,xlim=c(-0.10,0.10),ylim=c(0,200),xlab="Difference of Means")
  abline(v=obs,col="red")
  pval<-round(mean(abs(res)>=abs(obs)),3)
  text(x=0.05,y=150,label=paste("D=",round(obs,3),"\n p=",pval,sep = ""),cex=0.8)
  
  
  ## Reset panel
  par(mfrow=c(1,1))
  
  # ** Outcomes as ordinal variables ----
  
  # Run regressions and store results in a list
  fit1<-lm(id ~ treat_binary, data=data) # strength of ethnic id
  fit2<-lm(neighbor_relig ~ treat_binary, data=data) #ok with neighbor of diff religion
  fit3<-lm(neighbor_eth ~ treat_binary, data=data) #ok with neighbor of diff religion 
  fit4<-lm(neighbor_foreign ~ treat_binary, data=data) #ok with neighbor of diff nationality
  fit5<-lm(concern ~ treat_binary, data=data) # concerned about violent extremism
  fit6<-lm(message_binary ~ treat_binary, data=data) # recorded message
  fit7<-lm(trust_other_relig ~ treat_binary, data=data) #trust of other religious groups
  fit8<-lm(trust_other_eth ~ treat_binary, data=data)# trust of other ethnicities
  fit9<-lm(violence_justified ~ treat_binary, data=data) # violence justified
  reg_list<-list(fit1, fit2, fit3, fit4, fit5, fit6, fit7, fit8, fit9) 
  
  # Take coefficients and values of 95 and 90 Confidence intervals
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3), confint(fit3, level=.90))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4), confint(fit4, level=.90))
  fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5), confint(fit5, level=.90))
  fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6), confint(fit6, level=.90))
  fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7), confint(fit7, level=.90))
  fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8), confint(fit8, level=.90))
  fit9.df<-cbind(as.data.frame(fit9$coefficients), confint(fit9), confint(fit9, level=.90))
  
  # add labels
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Strength ethnic\nvs national id"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Neighbor\ndiff religion"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Neighbor\ndiff ethnicity"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Neighbor\nforeign"
  fit5.df$names<-rownames(fit5.df)
  fit5.df$outcome<-"Concern\nextremism"
  fit6.df$names<-rownames(fit6.df)
  fit6.df$outcome<-"Record\nmessage"
  fit7.df$names<-rownames(fit7.df)
  fit7.df$outcome<-"Trust\nother religion"
  fit8.df$names<-rownames(fit8.df)
  fit8.df$outcome<-"Trust\nother ethnicity"
  fit9.df$names<-rownames(fit9.df)
  fit9.df$outcome<-"Violence\njustified"
  
  # label table
  colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit5.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit6.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit7.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit8.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit9.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  
  #put everything together
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df, fit5.df, fit6.df, fit7.df, fit8.df, fit9.df)
  all_reg %<>% filter(!names=="(Intercept)") #take out intercept
  
  # Create plot
  all_reg %<>% mutate(outcome=factor(outcome, levels=unique(outcome)))
  plot<-all_reg %>% ggplot(aes(x=outcome, y=coef))  + geom_hline(yintercept=0, linetype=2, color="red") + geom_point(size=2) + geom_linerange(aes(ymin=lower_95, ymax=upper_95), size=.5) + geom_linerange(aes(ymin=lower_90, ymax=upper_90), size=1) + scale_y_continuous(limits=c(-1,1)) + labs(x="Outcome", y="Treatment effect")
  
  # Save
  save_plot("figA1.png", plot, base_height = 4, base_width = 7)
  
# ** Non-pooled Treatments ----
  
  
  
  # Run regressions and store results in a list
  fit1<-lm(id ~ treat, data=data) # strength of ethnic id
  fit2<-lm(neighbor_relig ~ treat, data=data) #ok with neighbor of diff religion
  fit3<-lm(neighbor_eth ~ treat, data=data) #ok with neighbor of diff religion 
  fit4<-lm(neighbor_foreign ~ treat, data=data) #ok with neighbor of diff nationality
  fit5<-lm(concern ~ treat, data=data) # concerned about violent extremism
  fit6<-lm(message_binary ~ treat, data=data) # recorded message
  fit7<-lm(trust_other_relig ~ treat, data=data) #trust of other religious groups
  fit8<-lm(trust_other_eth ~ treat, data=data)# trust of other ethnicities
  fit9<-lm(violence_justified ~ treat, data=data) # violence justified
  reg_list<-list(fit1, fit2, fit3, fit4, fit5, fit6, fit7, fit8, fit9) 

  # Take coefficients and values of 95  Confidence intervals
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4))
  fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5))
  fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6))
  fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7))
  fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8))
  fit9.df<-cbind(as.data.frame(fit9$coefficients), confint(fit9))
  
  # Label variables
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Ethnic vs national id"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Neighbor diff religion"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Neighbor diff ethnicity"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Neighbor foreign"
  fit5.df$names<-rownames(fit5.df)
  fit5.df$outcome<-"Concern about extremism"
  fit6.df$names<-rownames(fit6.df)
  fit6.df$outcome<-"Record message"
  fit7.df$names<-rownames(fit7.df)
  fit7.df$outcome<-"Trust other religion"
  fit8.df$names<-rownames(fit8.df)
  fit8.df$outcome<-"Trust other ethnicity"
  fit9.df$names<-rownames(fit9.df)
  fit9.df$outcome<-"Violence\njustified"
  
  # Label columns
  colnames(fit1.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit5.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit6.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit7.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit8.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit9.df)<-c("coef", "lower", "upper",  "names", "outcome")
  
  # Combine all reg results
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df, fit5.df, fit6.df, fit7.df, fit8.df, fit9.df)
  # Take out intercept
  all_reg %<>% filter(!names=="(Intercept)") 
  
  # Select and label treatments
  all_reg %>% select(names) %>% unique
  all_reg %<>% mutate(names=factor(names, levels=unique(names), labels=c("Scripture and ID", "Scripture and No ID"," Superordinate and ID"," Superordinate and No ID")))
  
  # Make outcomes into factor
  all_reg %<>% mutate(outcome=factor(outcome, levels=unique(outcome)))
  #Plot
  plot<-all_reg %>% ggplot(aes(x=outcome, y=coef, color=names))  + geom_hline(yintercept=0) + geom_point(size=1, position = position_dodge(width=.5)) + geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position = position_dodge(width=.5)) + scale_y_continuous(limits=c(-1,1)) + labs(x="Outcome", y="Treatment effect") + coord_flip() + scale_colour_ordinal("Treatment", begin=.2) 
  
  # Save plot
  save_plot("figA2.png", plot, base_height = 4, base_width = 7)
  
  
  # ** Muslims Only ----
  # The structure of the code is the same as the one above so we don't repeat annotations
  
  fit1<-lm(id ~ treat_binary, data=data_muslims) # strength of ethnic id
  fit2<-lm(neighbor_relig ~ treat_binary, data=data_muslims) #ok with neighbor of diff religion
  fit3<-lm(neighbor_eth ~ treat_binary, data=data_muslims) #ok with neighbor of diff religion 
  fit4<-lm(neighbor_foreign ~ treat_binary, data=data_muslims) #ok with neighbor of diff nationality
  fit5<-lm(concern ~ treat_binary, data=data_muslims) # concerned about violent extremism
  fit6<-lm(message_binary ~ treat_binary, data=data_muslims) # recorded message
  fit7<-lm(trust_other_relig ~ treat_binary, data=data_muslims) #trust of other religious groups
  fit8<-lm(trust_other_eth ~ treat_binary, data=data_muslims)# trust of other ethnicities
  fit9<-lm(violence_justified ~ treat_binary, data=data_muslims) # violence justified

  reg_list<-list(fit1, fit2, fit3, fit4, fit5, fit6, fit7, fit8, fit9)
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3), confint(fit3, level=.90))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4), confint(fit4, level=.90))
  fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5), confint(fit5, level=.90))
  fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6), confint(fit6, level=.90))
  fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7), confint(fit7, level=.90))
  fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8), confint(fit8, level=.90))
  fit9.df<-cbind(as.data.frame(fit9$coefficients), confint(fit9), confint(fit9, level=.90))
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Strength ethnic\nvs national id"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Neighbor\ndiff religion"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Neighbor\ndiff ethnicity"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Neighbor\nforeign"
  fit5.df$names<-rownames(fit5.df)
  fit5.df$outcome<-"Concern\nextremism"
  fit6.df$names<-rownames(fit6.df)
  fit6.df$outcome<-"Record\nmessage"
  fit7.df$names<-rownames(fit7.df)
  fit7.df$outcome<-"Trust\nother religion"
  fit8.df$names<-rownames(fit8.df)
  fit8.df$outcome<-"Trust\nother ethnicity"
  fit9.df$names<-rownames(fit9.df)
  fit9.df$outcome<-"Violence\njustified"
  
  colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit5.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit6.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit7.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit8.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit9.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df, fit5.df, fit6.df, fit7.df, fit8.df, fit9.df)
  all_reg %<>% filter(!names=="(Intercept)") 
  

  all_reg %<>% mutate(outcome=factor(outcome, levels=unique(outcome)))
  plot<-all_reg %>% ggplot(aes(x=outcome, y=coef))  + geom_hline(yintercept=0, linetype=2, color="red") + geom_point(size=2) + geom_linerange(aes(ymin=lower_95, ymax=upper_95), size=.5) + geom_linerange(aes(ymin=lower_90, ymax=upper_90), size=1) + scale_y_continuous(limits=c(-1,1)) + labs(x="Outcome", y="Treatment effect (Muslims only)")
  
  
  # Save plot
  save_plot("figA3.png", plot, base_height = 4, base_width = 7)
  
  # ** Failed Attention Check ----
  # Redo Figure 2(a) Attitudes toward outgroups
  # Regress each outcome on the binary treatment indicator and store the results
  fit1<-lm(religion_binary ~ treat_binary, data=data_attention_fail)
  fit2<-lm(ethnic_binary ~ treat_binary, data=data_attention_fail)
  fit3<-lm(trust_rel_binary ~ treat_binary, data=data_attention_fail)
  fit4<-lm(trust_eth_binary ~ treat_binary, data=data_attention_fail)
  reg_list<-list(fit1, fit2, fit3, fit4)
  # Grab the coefficients, 95%, and 90% CI of each and store in a dataframe
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3), confint(fit3, level=.90))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4), confint(fit4, level=.90))
  
  # Create columns that label treatment and outcomes 
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Tolerance of neighbor\nfrom  different religion"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Tolerance of neighbor\nfrom  different ethnicity"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Trusts members of \nother religious groups"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Trusts members of \nother ethnic groups"
  
  #Create labels for the remaining columns of dataframe
  colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  
  #Bring together the four regressions
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df)
  all_reg %<>% filter(!names=="(Intercept)") #keep only treatment effect
  all_reg$order <- seq(from = nrow(all_reg), to = 1, by = -1) # create a counter to reverse order of presentation
  
  #Create figure and store it
  ggplot(all_reg,aes(x=reorder(outcome,order),y=coef))+
    coord_flip()+
    geom_hline(yintercept = 0, color = "gray80", lty = "dashed")+
    geom_point(position=position_dodge(width = .5), size = 2.5) +
    theme(panel.background = element_blank(),
          legend.title = element_blank(), 
          plot.title = element_text(size = 12),
          panel.border = element_rect(colour = "gray", fill=NA, size=.8),
          legend.position = "bottom",
          axis.title.x  = element_text(size=10),
          axis.text.x  = element_text(angle=0, vjust=.5, hjust = .5, size=8))+
    geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
                  position=position_dodge(width = 0.5),
                  width=0, lwd = 0.4)+
    geom_errorbar(aes(ymin=lower_90, ymax=upper_90), 
                  position=position_dodge(width = 0.5),
                  width=0, lwd = 1) +
    scale_y_continuous(limits=c(-0.15,0.15)) + 
    labs(x="", y="Average Treatment Effect",title="Figure 2(a) Attitudes toward Outgroups (including attention check failures)")
  
  # Save it
  ggsave("figA4.png", width = 8, height = 6)

  
  # Redo Figure 2(b) Attitudes about extremism
  # Repeat same steps as above for attitudes about extremism
  fit1<-lm(concern_binary ~ treat_binary, data=data_attention_fail)
  fit2<-lm(justified_binary ~ treat_binary, data=data_attention_fail)
  
  reg_list<-list(fit1, fit2)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
  
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Concerned about\nviolent extremism"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Believes violence\n can be justified"
  
  colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  
  all_reg<-bind_rows(fit1.df, fit2.df)
  all_reg %<>% filter(!names=="(Intercept)") 
  all_reg$order <- seq(from = nrow(all_reg), to = 1, by = -1)

  ggplot(all_reg,aes(x=reorder(outcome,order),y=coef))+
    coord_flip()+
    geom_hline(yintercept = 0, color = "gray80", lty = "dashed")+
    geom_point(position=position_dodge(width = .5), size = 2.5) +
    theme(panel.background = element_blank(),
          legend.title = element_blank(), 
          plot.title = element_text(size = 10),
          panel.border = element_rect(colour = "gray", fill=NA, size=.8),
          legend.position = "bottom",
          axis.title.x  = element_text(size=8),
          axis.text.x  = element_text(angle=0, vjust=.5, hjust = .5, size=6))+
    geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
                  position=position_dodge(width = 0.5),
                  width=0, lwd = 0.4)+
    geom_errorbar(aes(ymin=lower_90, ymax=upper_90), 
                  position=position_dodge(width = 0.5),
                  width=0, lwd = 1) +
    scale_y_continuous(limits=c(-0.15,0.15)) + 
    labs(x="", y="Average Treatment Effect",title="Figure 2(b): Attitudes about Extremism (including attention check failures)")
  
  # Save it
  ggsave("figA5.png", width = 8, height = 6)
  
  # Redo Figure 2(c) Civic Participation
  # Repeat same steps as above for civic participation
  
  fit1<-lm(message_binary ~ treat_binary, data=data_attention_fail) #any message
  fit2<-lm(msg_gov~ treat_binary, data=data_attention_fail) # Message to government
  fit3<-lm(msg_islamic ~ treat_binary, data=data_attention_fail) # Message to Islamic council
  
  
  reg_list<-list(fit1, fit2, fit3)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3), confint(fit3, level=.90))
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Sends any message"
  
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Sends message to\ngovernment"
  
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Sends message to\nIslamic council"
  
  
  colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  
  
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df)
  all_reg %<>% filter(!names=="(Intercept)") 
  all_reg$order <- seq(from = nrow(all_reg), to = 1, by = -1)
  
  
  fig2c<-ggplot(all_reg,aes(x=reorder(outcome,order),y=coef))+
    coord_flip()+
    geom_hline(yintercept = 0, color = "gray80", lty = "dashed")+
    geom_point(position=position_dodge(width = .5), size = 2.5) +
    theme(panel.background = element_blank(),
          legend.title = element_blank(), 
          plot.title = element_text(size = 10),
          panel.border = element_rect(colour = "gray", fill=NA, size=.8),
          legend.position = "bottom",
          axis.title.x  = element_text(size=8),
          axis.text.x  = element_text(angle=0, vjust=.5, hjust = .5, size=6))+
    geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
                  position=position_dodge(width = 0.5),
                  width=0, lwd = 0.4)+
    geom_errorbar(aes(ymin=lower_90, ymax=upper_90), 
                  position=position_dodge(width = 0.5),
                  width=0, lwd = 1) +
    scale_y_continuous(limits=c(-0.15,0.15)) + 
    labs(x="", y="Average Treatment Effect",title="Figure 2(c) Civic Participation (including attention check failures)")
  
  
  # Save it
  ggsave("figA6.png", width = 8, height = 6)
  
  
  
  
  
  # ** Older Respondents Only ----
  # Redo Figure 2(a) Attitudes toward outgroups
  # Regress each outcome on the binary treatment indicator and store the results
  fit1<-lm(religion_binary ~ treat_binary, data=data_14_plus)
  fit2<-lm(ethnic_binary ~ treat_binary, data=data_14_plus)
  fit3<-lm(trust_rel_binary ~ treat_binary, data=data_14_plus)
  fit4<-lm(trust_eth_binary ~ treat_binary, data=data_14_plus)
  reg_list<-list(fit1, fit2, fit3, fit4)
  # Grab the coefficients, 95%, and 90% CI of each and store in a dataframe
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3), confint(fit3, level=.90))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4), confint(fit4, level=.90))
  
  # Create columns that label treatment and outcomes 
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Tolerance of neighbor\nfrom  different religion"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Tolerance of neighbor\nfrom  different ethnicity"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Trusts members of \nother religious groups"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Trusts members of \nother ethnic groups"
  
  #Create labels for the remaining columns of dataframe
  colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  
  #Bring together the four regressions
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df)
  all_reg %<>% filter(!names=="(Intercept)") #keep only treatment effect
  all_reg$order <- seq(from = nrow(all_reg), to = 1, by = -1) # create a counter to reverse order of presentation
  
  #Create figure and store it
  ggplot(all_reg,aes(x=reorder(outcome,order),y=coef))+
    coord_flip()+
    geom_hline(yintercept = 0, color = "gray80", lty = "dashed")+
    geom_point(position=position_dodge(width = .5), size = 2.5) +
    theme(panel.background = element_blank(),
          legend.title = element_blank(), 
          plot.title = element_text(size = 12),
          panel.border = element_rect(colour = "gray", fill=NA, size=.8),
          legend.position = "bottom",
          axis.title.x  = element_text(size=10),
          axis.text.x  = element_text(angle=0, vjust=.5, hjust = .5, size=8))+
    geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
                  position=position_dodge(width = 0.5),
                  width=0, lwd = 0.4)+
    geom_errorbar(aes(ymin=lower_90, ymax=upper_90), 
                  position=position_dodge(width = 0.5),
                  width=0, lwd = 1) +
    scale_y_continuous(limits=c(-0.15,0.15)) + 
    labs(x="", y="Average Treatment Effect",title="Figure 2(a) Attitudes toward Outgroups  (14 and older)")
  
  # Save it
  ggsave("figA7.png", width = 8, height = 6)
  
  
  # Redo Figure 2(b) Attitudes about extremism
  # Repeat same steps as above for attitudes about extremism
  fit1<-lm(concern_binary ~ treat_binary, data=data_14_plus)
  fit2<-lm(justified_binary ~ treat_binary, data=data_14_plus)
  
  reg_list<-list(fit1, fit2)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
  
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Concerned about\nviolent extremism"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Believes violence\n can be justified"
  
  colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  
  all_reg<-bind_rows(fit1.df, fit2.df)
  all_reg %<>% filter(!names=="(Intercept)") 
  all_reg$order <- seq(from = nrow(all_reg), to = 1, by = -1)
  
  ggplot(all_reg,aes(x=reorder(outcome,order),y=coef))+
    coord_flip()+
    geom_hline(yintercept = 0, color = "gray80", lty = "dashed")+
    geom_point(position=position_dodge(width = .5), size = 2.5) +
    theme(panel.background = element_blank(),
          legend.title = element_blank(), 
          plot.title = element_text(size = 10),
          panel.border = element_rect(colour = "gray", fill=NA, size=.8),
          legend.position = "bottom",
          axis.title.x  = element_text(size=8),
          axis.text.x  = element_text(angle=0, vjust=.5, hjust = .5, size=6))+
    geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
                  position=position_dodge(width = 0.5),
                  width=0, lwd = 0.4)+
    geom_errorbar(aes(ymin=lower_90, ymax=upper_90), 
                  position=position_dodge(width = 0.5),
                  width=0, lwd = 1) +
    scale_y_continuous(limits=c(-0.15,0.15)) + 
    labs(x="", y="Average Treatment Effect",title="Figure 2(b): Attitudes about Extremism (14 and older)")
  
  # Save it
  ggsave("figA8.png", width = 8, height = 6)
  
  # Redo Figure 2(c) Civic Participation
  # Repeat same steps as above for civic participation
  
  fit1<-lm(message_binary ~ treat_binary, data=data_14_plus) #any message
  fit2<-lm(msg_gov~ treat_binary, data=data_14_plus) # Message to government
  fit3<-lm(msg_islamic ~ treat_binary, data=data_14_plus) # Message to Islamic council
  
  
  reg_list<-list(fit1, fit2, fit3)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1), confint(fit1, level=.90))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2), confint(fit2, level=.90))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3), confint(fit3, level=.90))
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Sends any message"
  
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Sends message to\ngovernment"
  
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Sends message to\nIslamic council"
  
  
  colnames(fit1.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower_95", "upper_95", "lower_90", "upper_90", "names", "outcome")
  
  
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df)
  all_reg %<>% filter(!names=="(Intercept)") 
  all_reg$order <- seq(from = nrow(all_reg), to = 1, by = -1)
  
  
  fig2c<-ggplot(all_reg,aes(x=reorder(outcome,order),y=coef))+
    coord_flip()+
    geom_hline(yintercept = 0, color = "gray80", lty = "dashed")+
    geom_point(position=position_dodge(width = .5), size = 2.5) +
    theme(panel.background = element_blank(),
          legend.title = element_blank(), 
          plot.title = element_text(size = 10),
          panel.border = element_rect(colour = "gray", fill=NA, size=.8),
          legend.position = "bottom",
          axis.title.x  = element_text(size=8),
          axis.text.x  = element_text(angle=0, vjust=.5, hjust = .5, size=6))+
    geom_errorbar(aes(ymin=lower_95, ymax=upper_95), 
                  position=position_dodge(width = 0.5),
                  width=0, lwd = 0.4)+
    geom_errorbar(aes(ymin=lower_90, ymax=upper_90), 
                  position=position_dodge(width = 0.5),
                  width=0, lwd = 1) +
    scale_y_continuous(limits=c(-0.15,0.15)) + 
    labs(x="", y="Average Treatment Effect",title="Figure 2(c) Civic Participation (14 and older)")
  
  
  # Save it
  ggsave("figA9.png", width = 8, height = 6)
  
  
# Appendix Multple Comparison Corrections ----
  
  
  # ** Full Sample ----
  # Run all the regressions and collect the uncorrected p-values
  fit<-lm(religion_binary~treat_binary, data=data)
  p1<-as.numeric(summary(fit)$coefficients[2,4])

  fit<-lm(ethnic_binary~treat_binary, data=data)
  p2<-as.numeric(summary(fit)$coefficients[2,4])

  fit<-lm(trust_rel_binary~treat_binary, data=data)
  p3<-as.numeric(summary(fit)$coefficients[2,4])

  fit<-lm(trust_eth_binary~treat_binary, data=data)
  p4<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(concern_binary~treat_binary, data=data)
  p5<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(justified_binary~treat_binary, data=data)
  p6<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(message_binary~treat_binary, data=data)
  p7<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(msg_gov~treat_binary, data=data)
  p8<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(msg_islamic~treat_binary, data=data)
  p9<-as.numeric(summary(fit)$coefficients[2,4])
  
  # Concatenate all p-values 
  p<-round(c(p1,p2,p3,p4,p5,p6,p7,p8,p9),3)
  
  # Run all corrections
  p_bonferroni<- p.adjust(p, "bonferroni")
  p_BH <- round(p.adjust(p, "BH"),3)
  p_holm <- p.adjust(p, "holm")
  # Combine columns with unadjusted and corrected p-vals
  p_all<-cbind(p,p_bonferroni,p_BH,p_holm)
  p_all

  # Labels for each row
  labels_rows<-c("Tolerance different religion", 
            "Tolerance different ethnicity",
            "Trusts  other religion",
            "Trusts  other ethnicity ",
            "Concerned extremism",
            "Violence justified",
            "Sends any message",
            "Sends message to government",
            "Sends message to Islamic council")
  rownames(p_all)<-labels_rows
  #Labels for columns
  labels_cols<-c("Original p-values","Bonferroni  p-values","BH  p-values","Holm  p-values")
  colnames(p_all)<-labels_cols
  
  # Create table in LaTeX (highlights p-value that is p<0.05 in unadjusted)
  p_all %>%
    kbl(booktabs = TRUE,caption = "Multiple comparison corrections, full sample") %>%
    kable_styling(latex_options = c("scale_down","hold_position")) %>%
    row_spec(4, color = 'black', background = 'lightgray') 

  # ** Mossi only sample ----
  # Same structure but with non-Mossi only sample
  fit<-lm(religion_binary~treat_binary, data=data_mossi)
  p1<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(ethnic_binary~treat_binary, data=data_mossi)
  p2<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(trust_rel_binary~treat_binary, data=data_mossi)
  p3<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(trust_eth_binary~treat_binary, data=data_mossi)
  p4<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(concern_binary~treat_binary, data=data_mossi)
  p5<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(justified_binary~treat_binary, data=data_mossi)
  p6<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(message_binary~treat_binary, data=data_mossi)
  p7<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(msg_gov~treat_binary, data=data_mossi)
  p8<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(msg_islamic~treat_binary, data=data_mossi)
  p9<-as.numeric(summary(fit)$coefficients[2,4])
  
  p<-round(c(p1,p2,p3,p4,p5,p6,p7,p8,p9),3)
  
  p_bonferroni<- p.adjust(p, "bonferroni")
  p_bonferroni
  
  p_BH <- round(p.adjust(p, "BH"),3)
  p_BH
  
  
  p_holm <- p.adjust(p, "holm")
  p_holm
  
  p_all<-cbind(p,p_bonferroni,p_BH,p_holm)
  p_all

  labels_rows<-c("Tolerance different religion", 
                 "Tolerance different ethnicity",
                 "Trusts  other religion",
                 "Trusts  other ethnicity ",
                 "Concerned extremism",
                 "Violence justified",
                 "Sends any message",
                 "Sends message to government",
                 "Sends message to Islamic council")
  labels_cols<-c("Original p-values","Bonferroni  p-values","BH  p-values","Holm  p-values")
  rownames(p_all)<-labels_rows
  colnames(p_all)<-labels_cols
  
  p_all %>%
    kbl(booktabs = TRUE,caption = "Multiple comparison corrections, Mossi only sample") %>%
    kable_styling(latex_options = c("scale_down","hold_position")) %>%
    row_spec(4, color = 'black', background = 'lightgray')
  
  
  
  # ** Non-mossi ----
  # Same structure but with non-Mossi only sample
  fit<-lm(religion_binary~treat_binary, data=data_nomossi)
  p1<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(ethnic_binary~treat_binary, data=data_nomossi)
  p2<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(trust_rel_binary~treat_binary, data=data_nomossi)
  p3<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(trust_eth_binary~treat_binary, data=data_nomossi)
  p4<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(concern_binary~treat_binary, data=data_nomossi)
  p5<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(justified_binary~treat_binary, data=data_nomossi)
  p6<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(message_binary~treat_binary, data=data_nomossi)
  p7<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(msg_gov~treat_binary, data=data_nomossi)
  p8<-as.numeric(summary(fit)$coefficients[2,4])
  
  fit<-lm(msg_islamic~treat_binary, data=data_nomossi)
  p9<-as.numeric(summary(fit)$coefficients[2,4])
  
  p<-round(c(p1,p2,p3,p4,p5,p6,p7,p8,p9),3)
  p_bonferroni<- p.adjust(p, "bonferroni")
  p_BH <- round(p.adjust(p, "BH"),3)
  p_holm <- p.adjust(p, "holm")
  p_all<-cbind(p,p_bonferroni,p_BH,p_holm)

  
  labels_rows<-c("Tolerance different religion", 
            "Tolerance different ethnicity",
            "Trusts  other religion",
            "Trusts  other ethnicity ",
            "Concerned extremism",
            "Violence justified",
            "Sends any message",
            "Sends message to government",
            "Sends message to Islamic council")
  labels_cols<-c("Original p-values","Bonferroni  p-values","BH  p-values","Holm  p-values")
  rownames(p_all)<-labels_rows
  colnames(p_all)<-labels_cols
  
  p_all %>%
    kbl(booktabs = TRUE,caption = "Multiple comparison corrections, non-Mossi  sample") %>%
    kable_styling(latex_options = c("scale_down","hold_position"))

# Pre-Analysis Plan Tests ----
  # ** Shared Identity Hypotheses ----
  # H1a: Messaging emphasizing shared identity will improve attitudes toward outgroups relative to control.  
  # Regress each outcome on pooled treatment of shared identity (not scripture)
  # Religion tolerance outcome
  reg<-lm(religion_binary~treat_pool,data=data[data$treat_pool!="Scripture",])%>% # Store regression results
    tidy()
  reg$term[2]<-c("Treatment") # Call the independent variable Treatment for clarity

  # Output LaTeX table of treatment effect
  kable(reg,
        caption = "Testing  H1a for attitudes toward religious outgroups.",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")
  
  # Ethncity tolerance outcome
  reg<-lm(ethnic_binary~treat_pool,data=data[data$treat_pool!="Scripture",])%>%
    tidy()
  reg$term[2]<-c("Treatment")
  
  kable(reg,
        caption = "Testing  H1a for attitudes toward ethnic outgroups.",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")
  
  # Foreigner tolerance outcome
  reg<-lm(foreign_binary~treat_pool,data=data[data$treat_pool!="Scripture",])%>%
    tidy()
  reg$term[2]<-c("Treatment")
  
  kable(reg,
        caption = "Testing  H1a for attitudes toward immigrant outgroups.",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")

  
  # H1b: Messaging emphasizing shared identity will reduce support for extremist groups.  
  # Same idea/structure but use support for extremism outcomevars
  reg<-lm(concern_binary~treat_pool,data=data[data$treat_pool!="Scripture",])%>%
    tidy()
  reg$term[2]<-c("Treatment")
  
  kable(reg,
        caption = "Testing  H1b for effect on concern for extremism.",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")
  
  
  reg<-lm(justified_binary~treat_pool,data=data[data$treat_pool!="Scripture",])%>%
    tidy()
  reg$term[2]<-c("Treatment")
  
  kable(reg,
        caption = "Testing  H1b for effect on belief violence is justified.",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")
  
  
  #H1c: Messaging emphasizing shared identity will increase political participation. 
  # Same idea/structure but use messaging outcome
  reg<-lm(message_binary~treat_pool,data=data[data$treat_pool!="Scripture",])%>%
    tidy()
  reg$term[2]<-c("Treatment")
  
  kable(reg,
        caption = "Testing  H1c using indicator for whether respondent sent message or not.",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")

  
  
  # ** Scripture hypotheses----
  # H2a: Messaging invoking Quranic scripture will improve attitudes toward outgroups relative to control.  
  # Regress each outcome on pooled treatment of shared identity (not superordinate id)
  # Otherwise, same coding structure as above so we don't repeat annotations
  reg<-lm(religion_binary~treat_pool,data=data[data$treat_pool!="Superordinate",])%>%
    tidy()
  reg$term[2]<-c("Treatment")
  
  kable(reg,
        caption = "Testing  H2a for attitudes toward religious outgroups.",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")
  
  reg<-lm(ethnic_binary~treat_pool,data=data[data$treat_pool!="Superordinate",])%>%
    tidy()
  reg$term[2]<-c("Treatment")
  
  kable(reg,
        caption = "Testing  H2a for attitudes toward ethnic outgroups.",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")
  
  reg<-lm(foreign_binary~treat_pool,data=data[data$treat_pool!="Superordinate",])%>%
    tidy()
  reg$term[2]<-c("Treatment")
  
  kable(reg,
        caption = "Testing  H2a for attitudes toward immigrant outgroups.",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")
  
  
  # H2b: Messaging invoking Quranic scripture will decrease support for extremist groups relative to control.  
  reg<-lm(concern_binary~treat_pool,data=data[data$treat_pool!="Superordinate",])%>%
    tidy()
  reg$term[2]<-c("Treatment")
  
  kable(reg,
        caption = "Testing  H2b for effect on concern for extremism.",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")
  
  
  reg<-lm(justified_binary~treat_pool,data=data[data$treat_pool!="Superordinate",])%>%
    tidy()
  reg$term[2]<-c("Treatment")
  
  kable(reg,
        caption = "Testing  H2b for effect on belief violence is justified.",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")
  
  
  #H2c: Messaging invoking Quranic scripture will increase political participation.  
  reg<-lm(message_binary~treat_pool,data=data[data$treat_pool!="Superordinate",])%>%
    tidy()
  reg$term[2]<-c("Treatment")
  
  kable(reg,
        caption = "Testing  H2c using indicator for whether respondent sent message or not.",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")


  # ** Hypotheses about messages ----
  #H4:  Messaging emphasizing shared identity will increase respondents’ willingness to identify with one’s national identity relative to messages invoking Quranic scripture.  
    
  # To test this, regress binary variable of strength of national identity on pooled treatment
  
  reg<-lm(id_binary~treat_pool,data=data)%>%
    tidy()
  
  kable(reg,
        caption = "Testing  H4 using indicator for whether respondent id nationally or not.",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")
  

  # H5: Messaging involving Quranic scripture will make respondents more likely to turn to religious authority figures for action on peace and tolerance-promotion.  
  # To test this, regress likelihood to message Islamic council on pooled treatment indicator 
  
  reg<-lm(msg_islamic~treat_pool,data=data[data$treat_pool!="Superordinate",])%>%
    tidy()
  reg$term[2]<-c("Treatment")
  
  kable(reg,
        caption = "Testing  H5 using indicator for whether respondent sent message to Islamic leader",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")
  
    
  #H6: Messaging emphasizing shared identity will make respondents more likely to turn to national governmental bodies for action on peace and tolerance-promotion.  
  # To test this, regress likelihood to message government on pooled treatment indicator 
  reg<-lm(msg_gov~treat_pool,data=data[data$treat_pool!="Scripture",])%>%
    tidy()
  reg$term[2]<-c("Treatment")
  
  kable(reg,
        caption = "Testing  H6 using indicator for whether respondent sent message to government",
        col.names = c("Treatment", "Coefficient", "SE", "t", "p"),
        digits = c(0, 2, 3, 2, 3),
        align = c("l", "r", "r", "r", "r"),
        booktabs=T
  )%>%
    kable_styling(latex_options = "HOLD_position")
  
  
  
  # ** Heterogenous Treatment Effects----
  
  # H8: Treatment effects for messages invoking scripture will be larger for respondents that are more religious. 
  # Dataset of only those who pray often
  data_pray<-data %>% filter(!is.na(pray_freq))  %>% filter(pray_often==1)
  
  # Same analysis from before
  fit1<-lm(id_binary ~ treat_pool , data=data_pray[!data_pray$treat_pool=="Superordinate",])
  fit2<-lm(religion_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",]) #religion
  fit3<-lm(ethnic_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",]) #ethnicity
  fit4<-lm(foreign_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",]) #foreigners
  fit5<-lm(concern_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",])
  fit6<-lm(message_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",])
  fit7<-lm(trust_rel_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",])
  fit8<-lm(trust_eth_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",])
  reg_list<-list(fit1, fit2, fit3, fit4, fit5, fit7, fit8)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4))
  fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5))
  fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6))
  fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7))
  fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8))
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Ethnic vs\nnational id"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Neighbor\ndiff religion"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Neighbor\ndiff ethnicity"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Neighbor\nforeign"
  fit5.df$names<-rownames(fit5.df)
  fit5.df$outcome<-"Concern about\nextremism"
  fit6.df$names<-rownames(fit6.df)
  fit6.df$outcome<-"Record\nmessage"
  fit7.df$names<-rownames(fit7.df)
  fit7.df$outcome<-"Trust other\nreligion"
  fit8.df$names<-rownames(fit8.df)
  fit8.df$outcome<-"Trust other\nethnicity"
  
  colnames(fit1.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit5.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit6.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit7.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit8.df)<-c("coef", "lower", "upper", "names", "outcome")
  
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df, fit5.df, fit6.df, fit7.df, fit8.df)
  
  all_reg %<>% filter(!names=="(Intercept)") 
  
  all_reg %<>% mutate(outcome=factor(outcome, levels=unique(outcome)))
  
  # Store all results 
  effects_pray<-all_reg %>% ggplot(aes(x=outcome, y=coef)) + geom_hline(yintercept=0, linetype=2) + geom_point(size=1, position = position_dodge(width=.5)) + geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position = position_dodge(width=.5)) + scale_y_continuous(limits=c(-.5,.5)) + labs(x="", y="", title="Respondents pray frequently") + coord_flip() + scale_colour_ordinal("Treatment") + theme(legend.position = "none")

  # Do not pray often
  data_pray<-data %>% filter(!is.na(pray_freq))  %>% filter(pray_often==0)
  
  fit1<-lm(id_binary ~ treat_pool , data=data_pray[!data_pray$treat_pool=="Superordinate",])
  fit2<-lm(religion_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",]) #religion
  fit3<-lm(ethnic_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",]) #ethnicity
  fit4<-lm(foreign_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",]) #foreigners
  fit5<-lm(concern_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",])
  fit6<-lm(message_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",])
  fit7<-lm(trust_rel_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",])
  fit8<-lm(trust_eth_binary ~ treat_pool, data=data_pray[!data_pray$treat_pool=="Superordinate",])
  
  reg_list<-list(fit1, fit2, fit3, fit4, fit5, fit7, fit8)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4))
  fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5))
  fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6))
  fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7))
  fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8))
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Ethnic vs\nnational id"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Neighbor\ndiff religion"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Neighbor\ndiff ethnicity"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Neighbor\nforeign"
  fit5.df$names<-rownames(fit5.df)
  fit5.df$outcome<-"Concern about\nextremism"
  fit6.df$names<-rownames(fit6.df)
  fit6.df$outcome<-"Record\nmessage"
  fit7.df$names<-rownames(fit7.df)
  fit7.df$outcome<-"Trust other\nreligion"
  fit8.df$names<-rownames(fit8.df)
  fit8.df$outcome<-"Trust other\nethnicity"
  
  colnames(fit1.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit5.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit6.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit7.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit8.df)<-c("coef", "lower", "upper", "names", "outcome")
  
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df, fit5.df, fit6.df, fit7.df, fit8.df)
  all_reg %<>% filter(!names=="(Intercept)") 
  all_reg %<>% mutate(outcome=factor(outcome, levels=unique(outcome)))
  
  effects_no_pray<-all_reg %>% ggplot(aes(x=outcome, y=coef)) + geom_hline(yintercept=0, linetype=2) + geom_point(size=1, position = position_dodge(width=.5)) + geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position = position_dodge(width=.5)) + scale_y_continuous(limits=c(-.5,.5)) + labs(x="", y="", title="Respondents pray infrequently") + coord_flip() + scale_colour_ordinal("Treatment") + theme(axis.text.y = element_blank(), legend.position = "right")

  # Combine two plots 
  grid.arrange(effects_pray, effects_no_pray, nrow=1, ncol=2, widths = 5:4 )  

  # H9: Treatment effects for messages invoking super-ordinate identity will be larger for respondents that have had more contact with members of different out-groups.
  data_out<-data %>% filter(out_exposure==1)
  data_in<-data %>% filter(out_exposure==0)
  
  # Those that have contact with out group
  fit1<-lm(id_binary ~ treat_binary, data=data_out)
  fit2<-lm(religion_binary ~ treat_binary, data=data_out) #religion
  fit3<-lm(ethnic_binary ~ treat_binary, data=data_out) #ethnicity
  fit4<-lm(foreign_binary ~ treat_binary, data=data_out) #foreigners
  fit5<-lm(concern_binary ~ treat_binary, data=data_out)
  fit6<-lm(message_binary ~ treat_binary, data=data_out)
  fit7<-lm(trust_rel_binary ~ treat_binary, data=data_out)
  fit8<-lm(trust_eth_binary ~ treat_binary, data=data_out)
  reg_list<-list(fit1, fit2, fit3, fit4, fit5, fit7, fit8)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4))
  fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5))
  fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6))
  fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7))
  fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8))
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Ethnic vs\nnational id"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Neighbor\ndiff religion"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Neighbor\ndiff ethnicity"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Neighbor\nforeign"
  fit5.df$names<-rownames(fit5.df)
  fit5.df$outcome<-"Concern about\nextremism"
  fit6.df$names<-rownames(fit6.df)
  fit6.df$outcome<-"Record\nmessage"
  fit7.df$names<-rownames(fit7.df)
  fit7.df$outcome<-"Trust other\nreligion"
  fit8.df$names<-rownames(fit8.df)
  fit8.df$outcome<-"Trust other\nethnicity"
  
  colnames(fit1.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit5.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit6.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit7.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit8.df)<-c("coef", "lower", "upper", "names", "outcome")
  
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df, fit5.df, fit6.df, fit7.df, fit8.df)
  
  all_reg %<>% filter(!names=="(Intercept)") 
  
  all_reg %<>% mutate(outcome=factor(outcome, levels=unique(outcome)))
  
  het_out<-all_reg %>% ggplot(aes(x=outcome, y=coef))  + geom_hline(yintercept=0, linetype=2) + geom_point(size=1, position = position_dodge(width=.5)) + geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position = position_dodge(width=.5)) + scale_y_continuous(limits=c(-.5,.5)) + labs(x="", y="", title="Out-groups in neighborhood") + coord_flip() + scale_colour_ordinal("Treatment", begin=.2) + theme(legend.position = "none")
  
  # Those that have contact with in group only
  fit1<-lm(id_binary ~ treat_binary, data=data_in)
  fit2<-lm(religion_binary ~ treat_binary, data=data_in) #religion
  fit3<-lm(ethnic_binary ~ treat_binary, data=data_in) #ethnicity
  fit4<-lm(foreign_binary ~ treat_binary, data=data_in) #foreigners
  fit5<-lm(concern_binary ~ treat_binary, data=data_in)
  fit6<-lm(message_binary ~ treat_binary, data=data_in)
  fit7<-lm(trust_rel_binary ~ treat_binary, data=data_in)
  fit8<-lm(trust_eth_binary ~ treat_binary, data=data_in)
  reg_list<-list(fit1, fit2, fit3, fit4, fit5, fit7, fit8)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4))
  fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5))
  fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6))
  fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7))
  fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8))
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Ethnic vs\nnational id"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Neighbor\ndiff religion"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Neighbor\ndiff ethnicity"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Neighbor\nforeign"
  fit5.df$names<-rownames(fit5.df)
  fit5.df$outcome<-"Concern about\nextremism"
  fit6.df$names<-rownames(fit6.df)
  fit6.df$outcome<-"Record\nmessage"
  fit7.df$names<-rownames(fit7.df)
  fit7.df$outcome<-"Trust other\nreligion"
  fit8.df$names<-rownames(fit8.df)
  fit8.df$outcome<-"Trust other\nethnicity"
  
  colnames(fit1.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit5.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit6.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit7.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit8.df)<-c("coef", "lower", "upper", "names", "outcome")
  
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df, fit5.df, fit6.df, fit7.df, fit8.df)
  
  all_reg %<>% filter(!names=="(Intercept)") 
  
  all_reg %<>% mutate(outcome=factor(outcome, levels=unique(outcome)))
  
  het_in<-all_reg %>% ggplot(aes(x=outcome, y=coef))  + geom_hline(yintercept=0, linetype=2) + geom_point(size=1, position = position_dodge(width=.5)) + geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position = position_dodge(width=.5)) + scale_y_continuous(limits=c(-.5,.5)) + labs(x="", y="", title="Neighborhood only in-groups") + coord_flip() + scale_colour_ordinal("Treatment", begin=.2) + theme(axis.text.y = element_blank())
  
  grid.arrange(het_out, het_in, nrow=1, ncol=2, widths = 5:4)  

  # H10: Treatment effects for all messages will be larger for adolescents that have been previously exposed to violence.

  data_violence<-data %>% filter(violence_exp==1)
  fit1<-lm(id_binary ~ treat_binary, data=data_violence)
  fit2<-lm(religion_binary ~ treat_binary, data=data_violence) #religion
  fit3<-lm(ethnic_binary ~ treat_binary, data=data_violence) #ethnicity
  fit4<-lm(foreign_binary ~ treat_binary, data=data_violence) #foreigners
  fit5<-lm(concern_binary ~ treat_binary, data=data_violence)
  fit6<-lm(message_binary ~ treat_binary, data=data_violence)
  fit7<-lm(trust_rel_binary ~ treat_binary, data=data_violence)
  fit8<-lm(trust_eth_binary ~ treat_binary, data=data_violence)
  reg_list<-list(fit1, fit2, fit3, fit4, fit5, fit7, fit8)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4))
  fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5))
  fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6))
  fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7))
  fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8))
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Ethnic vs\nnational id"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Neighbor\ndiff religion"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Neighbor\ndiff ethnicity"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Neighbor\nforeign"
  fit5.df$names<-rownames(fit5.df)
  fit5.df$outcome<-"Concern about\nextremism"
  fit6.df$names<-rownames(fit6.df)
  fit6.df$outcome<-"Record\nmessage"
  fit7.df$names<-rownames(fit7.df)
  fit7.df$outcome<-"Trust other\nreligion"
  fit8.df$names<-rownames(fit8.df)
  fit8.df$outcome<-"Trust other\nethnicity"
  
  colnames(fit1.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit5.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit6.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit7.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit8.df)<-c("coef", "lower", "upper", "names", "outcome")
  
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df, fit5.df, fit6.df, fit7.df, fit8.df)
  
  all_reg %<>% filter(!names=="(Intercept)") 
  
  all_reg %<>% mutate(outcome=factor(outcome, levels=unique(outcome)))
  
  effects_violence<-all_reg %>% ggplot(aes(x=outcome, y=coef)) + geom_hline(yintercept=0, linetype=2) + geom_point(size=1, position = position_dodge(width=.5)) + geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position = position_dodge(width=.5)) + scale_y_continuous(limits=c(-.5,.5)) + labs(title="Exposed to violence", x="", y="") + coord_flip() + scale_colour_ordinal("Treatment") + theme(legend.position = "none") 
  
  
  data_no_violence<-data %>% filter(violence_exp==0)
  fit1<-lm(id_binary ~ treat_binary, data=data_no_violence)
  fit2<-lm(religion_binary ~ treat_binary, data=data_no_violence) #religion
  fit3<-lm(ethnic_binary ~ treat_binary, data=data_no_violence) #ethnicity
  fit4<-lm(foreign_binary ~ treat_binary, data=data_no_violence) #foreigners
  fit5<-lm(concern_binary ~ treat_binary, data=data_no_violence)
  fit6<-lm(message_binary ~ treat_binary, data=data_no_violence)
  fit7<-lm(trust_rel_binary ~ treat_binary, data=data_no_violence)
  fit8<-lm(trust_eth_binary ~ treat_binary, data=data_no_violence)
  reg_list<-list(fit1, fit2, fit3, fit4, fit5, fit7, fit8)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4))
  fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5))
  fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6))
  fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7))
  fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8))
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Ethnic vs\nnational id"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Neighbor\ndiff religion"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Neighbor\ndiff ethnicity"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Neighbor\nforeign"
  fit5.df$names<-rownames(fit5.df)
  fit5.df$outcome<-"Concern about\nextremism"
  fit6.df$names<-rownames(fit6.df)
  fit6.df$outcome<-"Record\nmessage"
  fit7.df$names<-rownames(fit7.df)
  fit7.df$outcome<-"Trust other\nreligion"
  fit8.df$names<-rownames(fit8.df)
  fit8.df$outcome<-"Trust other\nethnicity"
  
  colnames(fit1.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit5.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit6.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit7.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit8.df)<-c("coef", "lower", "upper", "names", "outcome")
  
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df, fit5.df, fit6.df, fit7.df, fit8.df)
  
  all_reg %<>% filter(!names=="(Intercept)") 
  
  all_reg %<>% mutate(outcome=factor(outcome, levels=unique(outcome)))
  
  effects_no_violence<-all_reg %>% ggplot(aes(x=outcome, y=coef)) + geom_hline(yintercept=0, linetype=2) + geom_point(size=1, position = position_dodge(width=.5)) + geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position = position_dodge(width=.5)) + scale_y_continuous(limits=c(-.5,.5)) + labs(title="No violence", x="", y="") + coord_flip() + scale_colour_ordinal("Treatment") + theme(axis.text.y = element_blank())
  
  grid.arrange(effects_violence, effects_no_violence, nrow=1, ncol=2, widths = 5:4)  


  # H11: Treatment effects for all religious messages will be larger for younger respondents.
  
  data_young<-data %>% filter(age<15)
  data_older<-data %>% filter(age>14)
  
  fit1<-lm(id_binary ~ treat_binary, data=data_young)
  fit2<-lm(religion_binary ~ treat_binary, data=data_young) #religion
  fit3<-lm(ethnic_binary ~ treat_binary, data=data_young) #ethnicity
  fit4<-lm(foreign_binary ~ treat_binary, data=data_young) #foreigners
  fit5<-lm(concern_binary ~ treat_binary, data=data_young)
  fit6<-lm(message_binary ~ treat_binary, data=data_young)
  fit7<-lm(trust_rel_binary ~ treat_binary, data=data_young)
  fit8<-lm(trust_eth_binary ~ treat_binary, data=data_young)
  reg_list<-list(fit1, fit2, fit3, fit4, fit5, fit7, fit8)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4))
  fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5))
  fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6))
  fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7))
  fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8))
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Ethnic vs\nnational id"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Neighbor\ndiff religion"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Neighbor\ndiff ethnicity"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Neighbor\nforeign"
  fit5.df$names<-rownames(fit5.df)
  fit5.df$outcome<-"Concern about\nextremism"
  fit6.df$names<-rownames(fit6.df)
  fit6.df$outcome<-"Record\nmessage"
  fit7.df$names<-rownames(fit7.df)
  fit7.df$outcome<-"Trust other\nreligion"
  fit8.df$names<-rownames(fit8.df)
  fit8.df$outcome<-"Trust other\nethnicity"
  
  colnames(fit1.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit5.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit6.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit7.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit8.df)<-c("coef", "lower", "upper", "names", "outcome")
  
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df, fit5.df, fit6.df, fit7.df, fit8.df)
  
  all_reg %<>% filter(!names=="(Intercept)") 
  
  all_reg %<>% mutate(outcome=factor(outcome, levels=unique(outcome)))
  
  het_younger<-all_reg %>% ggplot(aes(x=outcome, y=coef))  + geom_hline(yintercept=0, linetype=2) + geom_point(size=1, position = position_dodge(width=.5)) + geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position = position_dodge(width=.5)) + scale_y_continuous(limits=c(-.5,.5)) + labs(x="", y="", title="Age 10-14") + coord_flip() + scale_colour_ordinal("Treatment", begin=.2) + theme(legend.position = "none")
  
  
  fit1<-lm(id_binary ~ treat_binary, data=data_older)
  fit2<-lm(religion_binary ~ treat_binary, data=data_older) #religion
  fit3<-lm(ethnic_binary ~ treat_binary, data=data_older) #ethnicity
  fit4<-lm(foreign_binary ~ treat_binary, data=data_older) #foreigners
  fit5<-lm(concern_binary ~ treat_binary, data=data_older)
  fit6<-lm(message_binary ~ treat_binary, data=data_older)
  fit7<-lm(trust_rel_binary ~ treat_binary, data=data_older)
  fit8<-lm(trust_eth_binary ~ treat_binary, data=data_older)
  reg_list<-list(fit1, fit2, fit3, fit4, fit5, fit7, fit8)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4))
  fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5))
  fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6))
  fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7))
  fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8))
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Ethnic vs\nnational id"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Neighbor\ndiff religion"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Neighbor\ndiff ethnicity"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Neighbor\nforeign"
  fit5.df$names<-rownames(fit5.df)
  fit5.df$outcome<-"Concern about\nextremism"
  fit6.df$names<-rownames(fit6.df)
  fit6.df$outcome<-"Record\nmessage"
  fit7.df$names<-rownames(fit7.df)
  fit7.df$outcome<-"Trust other\nreligion"
  fit8.df$names<-rownames(fit8.df)
  fit8.df$outcome<-"Trust other\nethnicity"
  
  colnames(fit1.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit5.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit6.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit7.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit8.df)<-c("coef", "lower", "upper", "names", "outcome")
  
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df, fit5.df, fit6.df, fit7.df, fit8.df)
  
  all_reg %<>% filter(!names=="(Intercept)") 
  
  
  all_reg %<>% mutate(outcome=factor(outcome, levels=unique(outcome)))
  
  het_older<-all_reg %>% ggplot(aes(x=outcome, y=coef))  + geom_hline(yintercept=0, linetype=2) + geom_point(size=1, position = position_dodge(width=.5)) + geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position = position_dodge(width=.5)) + scale_y_continuous(limits=c(-.5,.5)) + labs(x="", y="", title="Age 15-18") + coord_flip() + scale_colour_ordinal("Treatment", begin=.2) + theme(axis.text.y = element_blank())
  
  grid.arrange(het_younger, het_older, nrow=1, ncol=2, widths = 5:4) 

# H12: Treatment effects for all messages from the imam will be larger for respondents that self-identify as Muslim compared to respondents that do not self-identify as Muslim.

  
  fit1<-lm(id_binary ~ treat_binary, data=data_muslims)
  fit2<-lm(religion_binary ~ treat_binary, data=data_muslims) #religion
  fit3<-lm(ethnic_binary ~ treat_binary, data=data_muslims) #ethnicity
  fit4<-lm(foreign_binary ~ treat_binary, data=data_muslims) #foreigners
  fit5<-lm(concern_binary ~ treat_binary, data=data_muslims)
  fit6<-lm(message_binary ~ treat_binary, data=data_muslims)
  fit7<-lm(trust_rel_binary ~ treat_binary, data=data_muslims)
  fit8<-lm(trust_eth_binary ~ treat_binary, data=data_muslims)
  reg_list<-list(fit1, fit2, fit3, fit4, fit5, fit7, fit8)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4))
  fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5))
  fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6))
  fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7))
  fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8))
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Ethnic vs\nnational id"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Neighbor\ndiff religion"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Neighbor\ndiff ethnicity"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Neighbor\nforeign"
  fit5.df$names<-rownames(fit5.df)
  fit5.df$outcome<-"Concern about\nextremism"
  fit6.df$names<-rownames(fit6.df)
  fit6.df$outcome<-"Record\nmessage"
  fit7.df$names<-rownames(fit7.df)
  fit7.df$outcome<-"Trust other\nreligion"
  fit8.df$names<-rownames(fit8.df)
  fit8.df$outcome<-"Trust other\nethnicity"
  
  colnames(fit1.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit5.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit6.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit7.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit8.df)<-c("coef", "lower", "upper", "names", "outcome")
  
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df, fit5.df, fit6.df, fit7.df, fit8.df)
  
  all_reg %<>% filter(!names=="(Intercept)") 
  
  all_reg %<>% mutate(outcome=factor(outcome, levels=unique(outcome)))
  
  effects_all_muslims<-all_reg %>% ggplot(aes(x=outcome, y=coef)) + geom_hline(yintercept=0, linetype=2) + geom_point(size=1, position = position_dodge(width=.5)) + geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position = position_dodge(width=.5)) + scale_y_continuous(limits=c(-.5,.5)) + labs(x="", y="", title="Muslim respondents") + coord_flip() + scale_colour_ordinal("Treatment") + theme(legend.position = "none")
  
  fit1<-lm(id_binary ~ treat_binary, data=data_nonmuslims)
  fit2<-lm(religion_binary ~ treat_binary, data=data_nonmuslims) #religion
  fit3<-lm(ethnic_binary ~ treat_binary, data=data_nonmuslims) #ethnicity
  fit4<-lm(foreign_binary ~ treat_binary, data=data_nonmuslims) #foreigners
  fit5<-lm(concern_binary ~ treat_binary, data=data_nonmuslims)
  fit6<-lm(message_binary ~ treat_binary, data=data_nonmuslims)
  fit7<-lm(trust_rel_binary ~ treat_binary, data=data_nonmuslims)
  fit8<-lm(trust_eth_binary ~ treat_binary, data=data_nonmuslims)
  reg_list<-list(fit1, fit2, fit3, fit4, fit5, fit7, fit8)
  
  fit1.df<-cbind(as.data.frame(fit1$coefficients), confint(fit1))
  fit2.df<-cbind(as.data.frame(fit2$coefficients), confint(fit2))
  fit3.df<-cbind(as.data.frame(fit3$coefficients), confint(fit3))
  fit4.df<-cbind(as.data.frame(fit4$coefficients), confint(fit4))
  fit5.df<-cbind(as.data.frame(fit5$coefficients), confint(fit5))
  fit6.df<-cbind(as.data.frame(fit6$coefficients), confint(fit6))
  fit7.df<-cbind(as.data.frame(fit7$coefficients), confint(fit7))
  fit8.df<-cbind(as.data.frame(fit8$coefficients), confint(fit8))
  
  fit1.df$names<-rownames(fit1.df)
  fit1.df$outcome<-"Ethnic vs\nnational id"
  fit2.df$names<-rownames(fit2.df)
  fit2.df$outcome<-"Neighbor\ndiff religion"
  fit3.df$names<-rownames(fit3.df)
  fit3.df$outcome<-"Neighbor\ndiff ethnicity"
  fit4.df$names<-rownames(fit4.df)
  fit4.df$outcome<-"Neighbor\nforeign"
  fit5.df$names<-rownames(fit5.df)
  fit5.df$outcome<-"Concern about\nextremism"
  fit6.df$names<-rownames(fit6.df)
  fit6.df$outcome<-"Record\nmessage"
  fit7.df$names<-rownames(fit7.df)
  fit7.df$outcome<-"Trust other\nreligion"
  fit8.df$names<-rownames(fit8.df)
  fit8.df$outcome<-"Trust other\nethnicity"
  
  colnames(fit1.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit2.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit3.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit4.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit5.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit6.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit7.df)<-c("coef", "lower", "upper", "names", "outcome")
  colnames(fit8.df)<-c("coef", "lower", "upper", "names", "outcome")
  
  all_reg<-bind_rows(fit1.df, fit2.df, fit3.df, fit4.df, fit5.df, fit6.df, fit7.df, fit8.df)
  
  all_reg %<>% filter(!names=="(Intercept)") 
  
  all_reg %<>% mutate(outcome=factor(outcome, levels=unique(outcome)))
  
  effects_all_nonmuslims<-all_reg %>% ggplot(aes(x=outcome, y=coef)) + geom_hline(yintercept=0, linetype=2) + geom_point(size=1, position = position_dodge(width=.5)) + geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position = position_dodge(width=.5)) + scale_y_continuous(limits=c(-.5,.5)) + labs(x="", y="", title="non-Muslim respondents") + coord_flip() + scale_colour_ordinal("Treatment")  + theme(axis.text.y = element_blank())

grid.arrange(effects_all_muslims, effects_all_nonmuslims, nrow=1, ncol=2, widths = 5:4)



